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ABSTRACT 

Feedback from supernovae is an essential aspect of galaxy formation. In order to improve 
subgrid models of feedback we perform a series of numerical experiments to investigate how 
supernova explosions shape the interstellar medium in a disk galaxy and power a galactic 
wind. We use the FLASH hydrodynamic code to model a simpUfied ISM, including grav- 
ity, hydrodynamics, radiative cooling above 10* K, and star formation that reproduces the 
Kennicutt-Schmidt relation. By simulating a small patch of the ISM in a tall box perpendic- 
ular to the disk, we obtain sub-parsec resolution allowing us to resolve individual supernova 
events. The hot interiors of supernova explosions combine into larger bubbles that sweep-up 
the initially hydrostatic ISM into a dense, warm cloudy medium, enveloped by a much hotter 
and tenuous medium, all phases in near pressure equilibrium. The unbound hot phase devel- 
ops into an outflow with wind speed increasing with distance as it accelerates from the disk. 
We follow the launch region of the galactic wind, where hot gas entrains and ablates warm 
ISM clouds leading to significantly increased mass loading of the flow, although we do not 
follow this material as it interacts with the galactic halo. 

We run a large grid of simulations in which we vary gas surface density, gas fraction, 
and star formation rate in order to investigate the dependencies of the mass loading, p = 
My,ind/M^,. In the cases with the most effective outflows we observe a /3 of 4, however in 
other cases we find /3 <C 1. We find that outflows are more efficient in disks with lower surface 
densities or gas fractions. A simple model in which the warm cloudy medium is the barrier 
that limits the expansion of the blast wave reproduces the scaling of outflow properties with 
disk parameters at high star formation rates. We extend the scaHng relations derived from an 
ISM patch to infer an effective mass loading for a galaxy with an exponential disk, finding that 
the mass loading depends on circular velocity as /3 cx V^" with a w 2.5 for a model which 
fits the Tully-Fisher relation. Such a scaling is often assumed in phenomenological models 
of galactic winds in order to reproduce the flat faint end slope of the mass function. Our 
normalisation is in approximate agreement with observed estimates of the mass loading for 
the Milky Way. The scahng we find sets the investigation of galaxy winds on a new footing, 
providing a physically motivated sub-grid description of winds that can be implemented in 
cosmological hydrodynamic simulations and phenomenological models. 
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1 INTRODUCTION 

Feedback is an essential aspect of galaxy formation models. It is in- 
voked to suppress the formation of large numbers of small galaxies 
(Rees & Ostriker 1977; White & Rees 1978; White & Frenk 1991). 
While photo heating can suppress star formation in the smallest ha- 
los, it cannot explain the low efficiency of SF in halos more massive 
than 10'' Mo (Efstathiou 1992; Okamoto et al. 2008). Feedback is 
also invoked to explain why such a small fraction of the baryons 
are in stars today (Fukugita et al. 1998; Balogh et al. 2001). An 
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efficient feedback implementation also appears essential for sim- 
ulations to produce realistic looking disk galaxies (Scannapieco 
et al. 2011; McCarthy et al. 2012). Observations of galactic winds 
at low (Heckman et al. 1990, 2000) and high z 3 redshift (Pet- 
tini et al. 2001) do show gas with a range of temperature and den- 
sities moving with large velocities of 100s of km s~^ with respect 
to the galaxy's stars, although the interpretation in terms of mass 
loss is complicated by the multi-phase nature of the wind (see e.g. 
Veilleux et al. 2005 for a recent review). Complimentary evidence 
for outflows comes from the high metal abundance detected in the 
IGM (Cowie et al. 1995), even at low densities (Schaye et al. 2003; 
Aguirre et al. 2004). Numerical simulations make it plausible that 
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galactic winds are responsible for this metal enrichment (Cen & 
Ostriker 1999; Aguirre et al. 2001 ; Theuns et al. 2002; Aguirre et al. 
2005; Oppenheimer & Dave 2006; Tescari et al. 2011), with low- 
mass galaxies dominating the enrichment of the bulk of the IGM 
(Booth et al. 2012). 

The sheer amount of energy released by supemovae (SNe) 
make the injection of energy into the interstellar medium (ISM) by 
SN explosions a prime candidate for driving galactic winds (Dekel 
& Silk 1986). However it is challenging to understand in detail how 
SNe regulate the transfer of mass and energy between the different 
phases of the ISM, as envisaged in the model of McKee & Ostriker 
(1977), and how and when this leads to the emergence of a galac- 
tic wind. Efstathiou (2000) and Silk (2001) extend the McKee & 
Ostriker (1977) model to examine how such interactions lead to 
self -regulation of star formation. They show that the properties of 
the galactic wind can be broadly understood once a temperature 
and density for the hot phase is found. This requires a model of 
evaporation of cold and warm clouds, yet without a more detailed 
understanding of the geometry and turbulence, we can go little fur- 
ther than steady spherically symmetric conduction models, which 
go back to Cowie & McKee (1977). Even if feedback is indeed due 
to SNe, it is not yet clear whether this is a consequence of their in- 
jection of hot gas, of turbulence (Mac Low & Klessen 2004; Scan- 
napieco & Briiggen 2010), of cosmic rays (Jubelgas et al. 2008), 
of the combined effects of magnetic fields, cosmic rays, and the 
galaxy's differential rotation (Kulpa-Dybel et al. 2011), or all of 
the above. 

Full hydrodynamic modelling of the interplay between the 
various components of the ISM in a Milky Way-like galaxy in 
a proper cosmological context is not yet currently possible due 
to the large range of scales involved, with density ranging from 
4 X 10~^^ g cm~* outside of halos to ~ 10"^° g cm"^ in cold 
clouds, temperatures from a few Kelvin inside star forming regions 
to ~ 10*K inside SN remnants, and time scales from a few thou- 
sand year for the propagation of a SN blast wave inside the ISM 
to ~ 10^° years for the age of the Galaxy. Excitingly, such full 
hydro-dynamical modelling begins to be possible for higher red- 
shift dwarfs (e.g Wise & Abel 2008), but for the moment models of 
larger galaxies at z ~ are limited to simulating a patch of galactic 
disk. In addition, we would also like to identify and understand the 
physics that is important in driving material from the galactic disk, 
and so it is desirable to have a series of numerical experiments. This 
is the approach we will follow in this paper. 

We begin by discussing constraints on galactic winds derived 
from current theoretical models of galaxy formation, and place our 
work in the context of comparable approaches. In section 3 we in- 
troduce the set-up of our own simulations. Briefly, we use a very 
simple model of the ISM which neglects the cold phase, and which 
is stirred by hot gas injected by SN explosions. Next we demon- 
strate that our sub-pc simulations have sufficient resolution to re- 
solve individual explosions, and illustrate the behaviour of both the 
ISM and of the wind for a reference model with properties chosen 
to be similar to that of the solar neighbourhood. In Section 5 we 
vary the properties of the simulated ISM (total and gas surface den- 
sities, star formation rate, cooling rate), and investigate if and when 
a wind is launched, and how its properties depend on that of the 
ISM. We obtain scaling relations of the wind to the ISM and ap- 
ply them in Section 6 to predict wind properties for a full galactic 
disk, and investigate how the wind properties depend on the galaxy 
properties. We summarise in Section 7. 



2 CONSTRAINTS ON GALACTIC WINDS 

2.1 Model requirements and observations 

We will assume that the baryon fraction in Milky-Way-sized halos, 
and halos of lower mass, falls significantly below the cosmological 
value, fb = Qb/^M due to the action of a galactic wind. Let the 
gaseous mass outflow rate from this wind be Mwind, and the star 
formation rate M*. A simple way to parameterise the efficiency of 
the SN-driven wind in removing baryons from the halo, is its mass 
loading, i.e. the ratio 
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where our P is equivalent to the /3 of Stringer et al. (2011). We 
introduce the hat in order to distinguish the average mass loading 

for a galaxy, B, from a local mass loading /3 at some point on the 
disk. If a galaxy exhausts its gas supply in star formation (and does 
not recycle wind material) then we will be left with a gas poor 
galaxy with baryon fraction reduced by a factor 1/{1 + fi). 

In order to infer the fraction of baryons ejected from galax- 
ies we can use the statistics of galaxies and dark matter halos. The 
number density of halos as a function of their mass can be approx- 
imated for masses below the exponential cut-off scale as a power 
law (Press & Schechter 1974; Reed et al. 2007), 
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Contrast this with the slope of the galaxy stellar mass function at 
low masses. 



(3) 



dlogM^, 

where observationally a is found to be in the range [—1.5, — 1], 
(see e.g. Blanton et al. 2003, 2005; Baldry et al. 2005, 2012; Li 
& White 2009). Naively identifying each dark matter halo with a 
galaxy of a given stellar mass (e.g. Guo et al. 2010) yields a galaxy 
mass to halo mass relation of M* oc M~°'^^^°''^^\ Identifying the 
stars as the main baryonic component implies a mass loading that 
scales with halo mass relatively steeply as (see also Stringer et al. 
2011) 

1 + ^ = /j, ^ oc m(i-«+«)/(i+") oc M-°-« , (4) 



where we substituted a faint end slope of a = —1.5 to derive the 
last exponent. Notably, this exponent — > oo as a — > — 1 and falls 
to zero as a — ^ —1.9, as such it is rather poorly constrained even 
by a well measured slope of the galaxy steUar mass function at low 
masses. One can infer not only that at low masses the mass loading 
/3 ^ 1 but also that it is strongly increasing towards lower-mass 
galaxies. 

Assume star formation results in the explosion of £ioo super- 
novae per 100 M0 of stars formed, each with energy Ssn, and that 
a fraction rjT gets converted into kinetic energy of an outflow. Ne- 
glecting other sources of energy then implies that 
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f]T £100 



10"'i org 



(5) 



where Vwind is the wind speed. If eioo, the thermalisation factor r)T 
and i?sN are all constants, then the product /3u,vind is also a con- 
stant. In this case large values of j3 imply lower wind speeds, and 
vice versa. If the mass-loading P indeed increases with decreasing 
halo mass, then of course eventually $ may become so large that the 
wind can no longer escape from the galaxy's potential well. Such 
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small halos may be subject to other destructive mechanisms, such 
as evaporation by re-ionization or obliteration by the explosions of 
the first stars. For massive halos, in order for the wind to escape it 
requires high wind speeds, implying low mass loading. The semi- 
analytical model of galaxy formation presented recently by Bower 
et al. (2012) imposes similar constraints on galactic winds to obtain 
fits to the faint-end of the galaxy mass function as inferred from our 
naive expectations: galactic winds need to have values of the mass 
loading /3 ~ 1 for Milky Way-like galaxies, with an indication 
that 3 increases even further towards lower masses. The best fitting 
models have /3 ~ 10 giving Vwind ~ 300 kms"'^. 

Numerical simulations of galaxy formation also try to imple- 
ment galactic winds with similar properties. Cosmological simula- 
tions such as Oppenheimer & Dave (2008) essentially implement 
the mass loading by hand by de-coupling the winds from the sur- 
rounding gas. More advanced techniques impose some constraints 
during the early stages of a burst of star formation when it is be- 
neath the simulation resolution, but later allow the gas distribution 
to evolve normally and the mass loading to emerge. Progress in 
this area has been made by simulations such as Dubois & Teyssier 
(2008) and Shen et al. (2010). Generally these include efficient 
feedback in an effort to produce a reasonable galaxy population, al- 
though they struggle to produce significant winds to remove enough 
baryons from Milky Way- sized galaxies. 

The OWLS simulations (Schaye et al. 2010) examined a va- 
riety of feedback prescriptions and models with efficient feedback 
in terms of a strong galactic winds fit a variety of properties of the 
galaxy population, including the TuUy-Fisher relation (McCarthy 
et al. 2012). However, in such models the properties of the winds 
are still part of the sub-grid modelling, i.e. the wind's properties are 
not computed but rather are simply imposed. This is required since 
the mass of gas entrained by a single supernova is a tiny fraction of 
the mass resolution element of the simulation (Creasey et al. 201 1). 

In order to directly simulate the generation of galactic winds 
requires a much higher resolution than can be reached in current 
cosmological simulations, as the sites of energy injection must be 
resolved (discussed further in section 3). In order to relax these con- 
straints, many simulators have either moved to high redshift (where 
the volumes are smaller), or modified the SNe in some way (such 
as aggregation of the energy injection). Examples of the former in- 
clude Mac Low & Ferrara (1999); Fujita et al. (2004); Wise & Abel 

(2008) and Powell et al. (2011), all of which struggle to produce 
mass loadings above unity except Wise & Abel (2008), who had 
massive Population in progenitors for their SNe. Examples of the 
latter include Dubois & Teyssier (2008) and Hopkins et al. (2011) 
with similar results, although Hopkins et al. (201 1) saw significant 
improvement by including the winds from massive stars. 

Despite having a different focus, there are also a number of 
studies of a high resolution SN driven ISM which have similar 
set-up to the current work, although they do not investigate the 
properties of their winds. loung & Mac Low (2006); Joung et al. 

(2009) ; Hill et al. (2012) and whilst this paper was being prepared 
Gent et al. (2012) have all modelled a SN driven ISM in a column 
through a galactic disk, driving a vertical wind. These studies in- 
vestigate the structure of the ISM, however their wind properties 
appear qualitatively similar to ours. On larger scales Strickland & 
Stevens (2000) and Cooper et al. (2008) have extended these to an 
approximation of the galaxy M82, although again the resolution 
restrictions severely limit the simulation time and SN energy injec- 
tion prescription. 

There are compelUng theoretical reasons to expect a high mass 
loading in galaxy winds, but are such winds seen in practise? The 



observational evidence for galactic outflows, at least in starhurst 
galaxies, is extremely strong (Heckman et al. 1990, 2000; Pettini 
et al. 2001; Martin 2005; Martin et al. 2002; Strickland & Heck- 
man 2009). Unfortunately it is notoriously difficult to constrain the 
wind properties from the data directly, partly because of uncertain 
metallicity and ionisation corrections needed to translate between 
the observed ion outflow and inferred total wind values, and partly 
because observing the wind in the spectrum of its galaxy does not 
provide spatial information of where the absorbing gas is located 
(Bouche et al. 2011, but see Wilman et al. 2005; Swinbank et al. 
2009 for a few cases of resolved studies of winds). The outflow- 
ing gas is Ukely multi-phase in nature, comphcating further the in- 
terpretation of the data. The picture for non starburst galaxies is 
even more complex, with Strickland & Heckman (2009) noting the 
lack of evidence for superwinds in such galaxies. As Chen et al. 
(2010) point out, however, the evidence for the high velocity out- 
flows come from blueshifted absorbers such as Na D that are trac- 
ing cooler material which is a fraction of the wind (or Mgll, for 
example Weiner et al. 2009 in the Deep2 galaxies). As far as it can 
be measured, the velocity of the outflow seems to be only weakly 
dependent on the SFR (Rupke et al. 2005). Probing the circimi- 
galactic medium around galaxies with a sight line to a background 
quasar allowed Bouche et al. (2011) to infer values of /3 = 2 — 3 
and wind speeds Wwind = 150 — 300 km s~^ for a set of L* galax- 
ies at redshift z ~ 0.1. They claim these wind speeds are in fact 
below the escape speed, and hence we may be observing a galactic 
fountain rather than a proper outflow. 

The picture of SNe as the driver of galactic winds also has con- 
sequences in terms of the metallicity of the galaxy. As SNe inject 
both metals and energy we expect and find a corresponding metal- 
licity deficit for low mass galaxies (Tremonti et al. 2004), interest 
in which goes back to Larson (1974). Both simple models (Peeples 
& Shankar 2011; Dayal et al. 2012) and simulations (Finlator & 
Dave 2008) show that galactic winds are an essential ingredient to 
obtain the observed mass-metallicity relation in galaxies. 

Summarising, observations provide strong evidence for the 
presence of galactic winds in star forming galaxies, but the param- 
eters of such winds are currently not tightly constrained. Models 
that make recourse to such winds to quench star formation require 
relatively high values of the wind's mass loading, /3 ~ 1 for MW- 
like galaxies, with 3 increasing for lower mass galaxies. But do 
SNe-driven winds indeed meet these requirements, and if they do, 
why? 

2.2 Resolving SNe in the ISM 

Ideally one would wish to probe the efficiency with which star for- 
mation can drive winds with simulations that self consistently in- 
cluded all the relevant physics, i.e. a full galaxy containing a star 
forming ISM, those stars subsequently redistributing their energy 
as type II SNe explosions, including outflows and cosmological 
infall. Unfortunately the range of scales involved in this problem 
makes such an approach currently computationally infeasible. To 
progress we must either truncate our resolution at some scale before 
we have fully resolved the physics, or to truncate our physics such 
that the available resolution becomes sufficient. The former route is 
one where we assume that we understand the physical processes to 
a certain degree and make our best effort at the calculation, forcing 
us to go deeply in to convergence studies. The latter is that of the 
numerical experiment where it is assumed that a certain amount of 
numerical calculation is possible and we make our best effort to in- 
clude the processes, requiring us to make full comparison with the 
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real Universe to test these assumptions (many simulations, are, of 
course, a mixture of these approaches). Our focus will be on the lat- 
ter case, that of the numerical experiment. We will also restrict our- 
selves to looking at the launch region of the galactic wind, where 
gas is expelled from the galaxy but not necessarily from the halo. 
This is consistent with what is needed to improve subgrid models 
in semi-analytic models and hydro simulations. 

The motivation for our choice of scale relates to the need to 
resolve individual SN blast waves as they sweep the ISM (as for 
example described by Cox 1972). Briefly, such explosions involve 
three distinct stages (e.g. Truelove & McKee 1999), beginning with 
the very early stage during which SNe ejecta expand almost freely 
into the ISM. As the ejecta sweep-up ISM preceded by a shock, 
eventually a reverse shock will run back into the ejecta, heating 
them to very high temperatures, signalling the start of the Sedov- 
Taylor stage (Sedov 1959; Taylor 1950). In both stages radiative 
cooling is negligible and consequently they can be described by 
similarity solutions, but the transition between them cannot. Finally 
at late times, the hot interior of shocked ejecta cools radiatively, and 
the swept-up shell of ISM and ejecta continue to 'snow-plough' fur- 
ther into the ISM, conserving momentum. Thornton et al. (1998) 
examine these last two states using a set of 1 dimensional simula- 
tions of the evolution of explosions in a uniform ISM, examining in 
detail the transition from the ST-phase to the snow-plough phase. 
They claim that radiative cooling is efficient enough that typically 
only 10 per cent of the initial blast energy is transferred to the ISM. 
Notably, the amount of gas heated by these explosions is not a lin- 
ear function of the SN energy, indeed it is sub-linear, and thus we 
may expect that aggregating the energy injection of many SNe into 
a single event will underestimate the amount of gas heated. 

We would in principle like to resolve the earliest phase of the 
explosions when ejecta dominate, but in this paper we restrict our- 
selves to initiate our SNe in the Sedov-Taylor phase. The transi- 
tion between ejecta-dominated and ST-phase occurs approximately 
when the shock has swept-up an amount of of ISM mass that is 
comparable to that originally ejected. In low density regions the 
size of the bubble where the transition happens may then be rela- 
tively large, and it would be worthwhile investigating whether this 
matters; we intend to do so in future work. Given this limitation, 
and for the typical ISM densities near the centre of the disk in our 
simulations, it then suffices to resolve scales of order of a few par- 
sees to fully capture the cooling of the swept-up shell of ISM (e.g. 
Cox 1972), and such a simulation will be able to resolve both the 
cooling and some part of the adiabatic phase of the remnant. 

As such the dependence of our question upon sub-parsec phe- 
nomena can be seen only in two key areas, raising the following 
questions 

(i) Star formation occurs on these scales, and thus controls the 
distribution (in time and space) of supernovae. Does this affect the 
properties of the galactic wind, for example because supernovae 
explode in high density environments and/or near to other super- 
novae? 

(ii) The medium that the SNe drive into contains structures on 
sub-parsec scales, for example cores of molecular clouds. Does this 
departure from a classical fluid affect the large scale wind? 

We will argue that the answers to both the above the questions may 
indeed be negative, motivating a set of simulations of a highly sim- 
plified ISM. Such a simulation would also improve our physical 
understanding of the role of the individual processes. 

On the first question we note that the progenitor of type II 
core collapse SNe are massive stars (Smartt 2009) with lifetimes 



~ 1 — 30 Myrs (Portinari et al. 1998), therefore the majority of 
SN energy associated with an instantaneous burst of stars with for 
example a Chabrier (2003) initial mass fimction will be released 
within ~ 30 Myrs. It is thought the birth cloud of such stars is 
likely destroyed before by the combined effects of stellar winds, 
proto-stellar jets and radiation (e.g. Matzner 2002), and there is ob- 
servational evidence for this (e.g. Lopez et al. 2011). Some clouds 
may form by turbulent compression when overrun by a spiral arm, 
and may disperse by the same flows that created them in the first 
place on a short time-scale (Dobbs 2008, see also Tasker & Tan 
2009). 

In any case, when the SNe explodes it will in general not do 
so inside its natal cloud. For this reason we assume that the SNe 
explode in typical environments in the disk plane of galaxies. Note 
however that the SNe may still be clustered rather than Poisson, a 
complication that we neglect. Typical giant molecular clouds have 
a velocity dispersion of ~ 4 km (Scoville & Sanders 1987; 
McCray & Kafatos 1987), which over 10 Myr results in a disper- 
sion of around 40 pc, which is a significant fraction of our box size 
and the typical distance between molecular clouds. 

The second question is deficate, and worthy of significant dis- 
cussion. We first note that we follow the nomenclature of Wolfire 
et al. (1995), where the T ~ lOOK phase of the ISM is called the 
cold neutral medium (hereafter CNM), the T ~ IO^'K phase as the 
warm neutral mediimi (hereafter WNM) and the T > lO^K phase 
as the hot ionised medium (HIM). The CNM exists in the form of 
dense clouds, occupying a very small fraction of the total volume 
but with a significant fraction of the total mass. These clouds are be- 
lieved to be in pressure equilibrium Spitzer (1956), with the WNM 
and HIM, thus making their energy budget (pressure x volume) 
also a small fraction of the ISM thermal energy. Their pressure sup- 
port is probably composed of a combination of magnetic, thermal 
and cosmic ray components. The proportions of energy in thermal, 
bulk and turbulent motions of the HIM and WNM are still not en- 
tirely known though there is consensus that much of the turbulence 
is supersonic (Elmegreen & Scalo 2004). A supersonic nature of 
turbulence in the ISM requires that the energy budget is dominated 
by inertial terms of the turbulent motions over the thermal and mag- 
netic terms in the WNM and HIM. 

Despite their small fraction of the energy budget, however, the 
cold phase can perform the role of a heat sink. Thermal energy from 
the warm and hot phases can be transported in to the cold phase 
via thermal conduction which can be dissipated via the molecular 
transitions of this cold gas (particularly CO, H2), metal lines (in 
particular CI*), and dust. The excited states of the molecules, how- 
ever, are rather long lived and whilst they are certainly important for 
star formation they may not significantly cool the WNM phase of 
the ISM due to its sparse nature (de Jong et al. 1980; Martin et al. 
1996). The molecules also play an important role as an absorber 
of photo-ionizing radiation, however we will ignore radiative driv- 
ing here. The simulations described in this paper simply neglect 
the cold phase, by tnmcating the cooling function below a value of 
To = 10'' K. If we were to include cooling below To we would 
have to include significantly more physics (magnetic fields, heat 
conduction, diffuse heating): here we want to investigate and un- 
derstand the simpler yet still complex case of a two-phase medium. 

We have also intentionally left out the physics of cosmic rays 
(see, e.g. Pfrommer et al. 2007) and magnetic fields (e.g. Bre- 
itschwerdt & de Avillez 2007) which may be important in provid- 
ing support against collapse, particularly in the CNM. Our goal is 
to imderstand the resultant ISM without these complications, be- 
fore discussing the implications of their addition. We would also 
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like to stress that although we have included gravity, we have not 
included ie/f -gravity (i.e our gravity is time-independent and only 
self-consistent for the initial set-up) which would be a poor approx- 
imation if we had included the dense, cold material of the CNM. 
Without the CNM gravity does not influence material on scales be- 
low the Jeans length of the WNM, equivalent to the scale height of 
the warm disk. 



3 SIMULATION SET-UP 

In the following section we will describe the simulations we have 
performed of supernova driven outflows from an idealised ISM. 
Our simulations model the ISM and halo of a disk galaxy in a tall 
column, with long (z) axis perpendicular to the galactic disk, and 
co-rotating with the disk material. We use outflow conditions at the 
top and bottom of the colunm, and periodic boundary conditions 
in X and y. We describe the initial conditions of the gas and the 
physical processes (gravity, cooling and supernova feedback) we 
have included, and detail their numerical implementation. Finally 
we describe some tests we have performed on the code and the 
parameters we chose to explore in our simulation set. 

The simulations use a modified version of the FLASH 3 code 
(Fryxell et al. 2000). FLASH 3 is a parallel, block structured, uni- 
form time-step, adaptive mesh refinement (AMR) code. Its second 
order (in space and time) scheme uses a piecewise-parabolic recon- 
struction in cells. Due to the extremely turbulent nature of the ISM 
in our simulations, we find that FLASH attempts to refine (i.e. to 
use the highest resolution allowed) almost everywhere within our 
simulation volume. Therefore we disable the AMR capability of 
FLASH and run it at a constant refinement level (albeit varied for 
our resolution studies). To mitigate the overhead of the guard-cell 
calculations we increase our block size to 32^ cells per block. 

For the gas physics we have assumed a monatomic ideal gas 
equation of state. 



(6) 



where u is the specific thermal energy and 7 = 5/,3 is the adi- 
abatic index. This deviates slightly from the physical equation of 
state which should include the ttansition in mean particle mass that 
occurs as the atomic hydrogen becomes ionised, but the impact of 
this simplification is small compared to the other uncertainties in 
this kind of simulation. 



3.1 Physical processes 

The simplified ISM discussed in Section 2.2 is shaped by three fun- 
damental processes: gravity, cooling and energy injection from su- 
pemovae, which dominate when we are only considering the WNM 
and HIM. We stress that our aim is to simplify the problem as much 
as possible so that we can extract the physical principles. In future 
works we will experiment with making some assumptions more 
realistic. Below we discuss the effects and implementation of all 
these processes. 



3.1.1 Gravity 

The gas in our simulations is initially in (vertical) hydrostatic equi- 
librium. In a disk galaxy the gravitational acceleration is induced 
by the gas and stars in the disk, baryons in the bulge and dark mat- 
ter (in the halo and possibly the disk, see e.g. Read et al. 2008). 
Despite these complications, when one moves to the (non-inertial) 



frame moving locally with the disk, the dominant effective poten- 
tial lies in the vertical direction, with a scale height of a few hun- 
dreds of parsecs. Since the shape of this profile is approximately 
in accordance with the gaseous one, we model the total gravity of 
all components (gas, stars, dark matter) as being in proportion to 
the gaseous component, with a multiplier of the inverse of the gas 
fraction, j- , to account for the steUar and dark matter components, 
i.e. the gravitational potential depends on the gas density through 
Poisson's equation as 



/gV> = 4nGp . 



(7) 



We also make a second assumption, namely that the gravita- 
tional profile of the disk is fixed in time, <j) = (jifpo]. This is as- 
sumed because the minimum temperature of our cooling function 
(discussed in Section 3.1.2 below) sets the Jeans length on the or- 
der of the scale of the disk height, so we do not expect smaUer 
self-gravitating clouds to appear in our simulations. In contrast, 
in the ISM of the Milky Way small self-gravitating clouds can 
form, because the ISM does cool to lower temperatures, however 
the physics of star formation is not the process we wish to address 
in this paper. 

Other terms we have neglected include those introduced by the 
Coriolis force across our simulation voltune, due to our non-inertial 

choice of frame. 



-2n A V , 



(8) 



where n is the angular velocity of the galaxy. Our simulations, 
however, will typically be of such short time scales and volumes 
that the Rossby number (the ratio of inertial to coriolis terms) is 
large. Nevertheless, more complete simulations would include this, 
along with the time dependent gravitational changes introduced by 
spiral density waves. Note that our simulations also neglect the ve- 
locity shear that is present in a differentiaUy rotating disc. 

5.7.2 Radiative cooling 

The cooling function A(r) of T ~ 10'' - lO'^K gas with solar 
abundances is primarily due to bound-bound and bound-free tian- 
sitions of ions, whereas above T — lO^K it is largely dominated by 
bremsstrahlung (Sutherland & Dopita 1993). Below T ~ 8000K 
there is a sharp decrease by several orders of magnitude, causing a 
build up of gas in the WNM. Cooling below 8000K is due to dust, 
metal transition lines such as CI*, and at very low temperatures, 
molecules. 

Whilst the imprint of small features in the cooling function 
should be observable in the ISM, it is really the cut-off at T ~ 
8000 K that controls the WNM, and as such we approximate the 
cooling function with a Heaviside function with a step at To = 
10* K, 



pu ■■ 



-An^ T^To 



0, 



T <To., 



(9) 



where we in addition assume pure hydrogen gas so that the number 
density n = p/rrip, and A = 10^^^ erg cm'* s^^ (although it is 
varied in a few of the simulations). We implement this very simple 
functional form so that we can expUcifly check the effect of the 
normalisation of the cooling function, and to make sure that any 
characteristic temperature of the gas is not due to features in A. 

The cooling function of Eq. (9) results in a cooling time for 
gas with T ^ To of 
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in 22 3 1 ) • (10) 

10~^-^ erg crri'^ y 

Since we have chosen a discontinuous function for our cooling, 
we implement a scheme in our code which prevents cooling below 
To (although the hydrodynamic forces can still achieve lower tem- 
peratures adiabatically). This largely prevents the overshoot errors 
resulting from an explicit solver in this kind of problem. 

To test the importance of the choice of cooling function, we 
also implemented the cooling function appropriate for cosmic gas 
with solar abundance pattern from Sutherland & Dopita (1993), 

AsD(r) = 5.3x10-"'' (rgi/^ O-S/^Tg-i/") erg erne's-' ,(11) 

where Tg = T/IO* K, with = 0.03 for low metallicity gas, 
and A = for T < lO"* K. All runs where this cooling function is 
used are marked 'SD' (see table 1). The minimum of this cooling 
function is at 5 x 10^/m K, where the cooling rate 



AsD.min = 1.30 X 10 erg CIl/ s ^ , 



(12) 

(ignoring the cut-off below 10* K). We show in Appendix A that 
the behaviour of the ISM in our simulations is surprisingly inde- 
pendent on the exact shape of the cooling function A(T), although 
it depends on the minimum value at high temperatures » 10* K. 



3.1.3 Energy injection by supemovae 

The Kennicutt-Schmidt (KS) relation connects observed surface 
density of star formation in a disk galaxy, E,, to its gas surface 
density Eg, 

w 2.5 X 10-*eJi* Mq yr"^ kpc"" , (13) 

(Kennicutt 1998), where Sgi = Eg/1 MqPc"". We also perform 
some simulations with an alternative formulation using a higher 
star formation rate, more commonly used in cosmological simula- 
tions, discussed in Appendix A. Notably this introduces an addi- 
tional dependence on the gas fraction of the disk, /g, that is absent 
from the KS relation. Our idealised model of a supernova explo- 
sion is the injection of 10''^ ergs (Cox 1972) of thermal energy in a 
small volume, implicitly assuming instantaneous thermalisation of 
the SN ejecta. The distribution in time of these is taken to be a Pois- 
son process (the Poisson process has the Markov property and so 
our SNe are independent) with a time independent rate computed 
from the initial parameters of the disk. For the local spatial distri- 
bution of SNe we assume the star formation rate to be proportional 
to the initial density, i.e. 



E[p* dV dt] = E* p!^_3. dV dt . 



(14) 



A consequence of this choice is that if the scale height of the gas 
profile evolves significantly the distribution of SNe will become 
inconsistent with the instantaneous mass profile. We discuss this 
further later. 

Given the star formation rate, the associated core-collapse 
SN rate is computed assuming the stellar initial mass function 
yields eioo SNe per 100 M© of star formation. For reference, for 
a Chabrier (2003) initial mass function with stars with masses e 
[0.1. 100] Mq, of which those with mass in the range [6, 100] M© 
undergo core collapse, eioo = 1-8. 

The final element of the SN prescription is the distribution 
of the injected energy over the computational grid. The choice of 



volume over which to spread the thermal energy of the supernovae 
is influenced by two considerations. If the volume is too large the 
remnant will evade the adiabatic expansion phase and immediately 
proceed to the radiative phase (Cox 1972; Creasey et al. 2011). If 
the volume is very small the code will require many extra time 
steps evaluating the initial stages of a Sedov-Taylor blast wave and 
will perform unnecessary computation^ . Following Cox (1972), the 
radius at which the blast wave cools and forms a dense shell is 
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(15) 



however to account for higher densities and the numerical spread- 
ing of shocks it is wise to resolve a fraction of this (Creasey et al. 
2011). 

Taking the above into consideration, for our simulations we 
spread the thermal energy of each SN over several cells given by 
the multivariate (3D) normal distribution of standard deviation 2 
pc, consistent with being smaller than the cooling radius of Eq. (15) 
for densities n < 77cm~^ (p < 1.3 x 10""" g cm~^). 



3.1.4 Time-stepping 

In addition to the numerical considerations described above, we 
also needed to make some adjustments to the time step calcula- 
tion in FLASH. The default time-stepping scheme in FLASH uses 
a Strang-split method (Strang 1968, an operator splitting method 
where the hydrodynamic update occurs in two half steps, with the 
order in which the Riemann-solver operates reversed from xyz to 
zyx between the first and the second half step). Source terms such 
as the injection of SN energy, are evaluated at the end of each half 
step, after the Riemann solver has been applied. This makes the im- 
plementation of the supernova energy injection problematic, as the 
thermal energy in a cell can increase by many orders magnitude fol- 
lowed by a hydrodynamic step before a new time step is calculated. 
The latter hydrodynamic step then almost inevitably violates the 
CFL condition and the Riemann solver fails to converge. We avoid 
this by making the timestep limiter for the supernova source terms 
predictive, i.e. we utilise the foreknowledge of the pre-computed 
SNe times to recognise when a supernova will occur before the end 
of the timestep given by the CFL condition and return a timestep 
of either up to just before the supernova, or of the predicted CFL 
timestep after the supernova has occurred, whichever is smaller. 

It is worth contrasting this with some other simulations of the 
ISM. In a series of simulation de Avillez & Breitschwerdt (2004, 
2005a,b) uses a set up similar to ours, with imposed gravity, cool- 
ing, SNe turbulence and magnetic fields in columns through disks 
of 1 X 1 X 10 kpc, although the focus is not on the mass loading. 
More recently the ERIS simulations (Powell et al. 201 1) simulated 
the ISM in a single high redshift dwarf galaxy. Cooper et al. (2008) 
perform a simulation of the central region of an M82-like starburst 
galaxy with gravity, cooling and energy injection due to supemovae 
(although this energy injection is continuous within a central vol- 
ume, rather than stochastic as in our simulations). 



^ To get some idea of the computational requirement of this, we recall that 
the velocity of a 3 dimensional Sedov blast wave evolves as « ~ t~^l^. 
Substituting this into the Courant-Friedrichs-Lewy (CFL) condition we see 
that the number of time steps required to reach a given radius is proportional 
to that radius. 
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3.1.5 Code tests 

A set up as complex as this requires some testing to confirm that 
the physical processes have been correctly implemented. As such 
we ran a number of simpler problems as well as the convergence 
tests in Appendix A. 

In order to test our hydrostatic set up we simulated the disk 
without supernovae for several dynamical times. Some sub-percent 
evolution in the gas occurred, almost certainly due to our evaluation 
of the analytic solution for the gravitational potential and density at 
the centres of cells producing some discretisation error. The im- 
plementation of the cooling function was tested largely in Creasey 
et al. (2011). We follow a similar approach where we made the 
cooling rate for each cell an output of our code which was com- 
pared with the instantaneous rate predicted from the temperature 
and density of each cell (again there were small differences due to 
the comparison of an instantaneous rate with the average from an 
implicit scheme). 

The implementation of the individual SN in our set-up is 
largely similar to that of the Sedov-Taylor blast wave solution im- 
plemented in FLASH as a standard test, and we compared it to the 
similarity solution. We calculate the location and times of SNe ex- 
plosions ahead of the simulation, and verify that the code indeed 
injects them correctly. 

We initially also performed these calculations using the GAD- 
GET simulation code (Springel 2005) that has been successfully 
applied to many cosmological simulations. Unfortunately the adap- 
tive time-stepping algorithm proved problematic for correctly fol- 
lowing the blast waves, and we noticed similar problems as recently 
highlighted by Durier & Dalla Vecchia (2012): particles may be on 
long time-steps in the cold ISM, and largely fail to properly ac- 
count for being shocked by the blast wave from a nearby particle. 
Durier & Dalla Vecchia (2012) addressed this problem with a time 
step propagation algorithm, however we did not have this nor the 
algorithm of Saitoh & Makino (2009) available and the alternative 
of a global timestep would have been far too computationally ex- 
pensive due to the large dynamic range in time steps required in the 
evolution of the blasts. As such we used the global adaptive time 
stepping algorithm of FLASH. 



3.2 Initial conditions 

Our initial setup is a tall box poking vertically through an idealised 
disk profile. We choose the long axis in the z-direction in order 
to capture a multiple of the gravitational scale height of the disk. 
The profile is a 1 -dimensional gravitationaUy bound isothermal one 
with gas surface density Eg. As discussed in section 3.1 we have 
excluded the effects of shear (due to the Coriolis force in the disk) 
and large scale motions which may drive some turbulence down to 
the small scales. The gas density is 

Ea 



(16) 



and the corresponding gravitational acceleration follows from 
Eq. (7), 



V$ = 27rGEg/g-Hanh (^^^ 



(17) 



Setting the gas temperature to To (which is also the base of the 
imposed cooling function) and assuming the gas to be initially in 
hydrostatic equilibrium, the scale height is 





Range of values explored 


Fiducial 
value 


Eg (M0 pc-2) 


2.5, 3.23,4.17,5.39, 6.96, 


11.61 




8.99, 11.61, 15, 30, 50, 150, 500 




/g 


0.01,0.015,0.022,0.033, 






0.050,0.1,0.2,0.5, 1.0 


0.1 




Eq. (13), (Al) 


Eq. (13) 


A (erg cm'^s-^) 


1,2,4, 8, 16x10-22, SD 


10-22 


Resolution (pc) 


0.78, 1.56,3.12, 6.25 





Table 1. Parameter variations in our simulation. Each simulation is ini- 
tialised with an isothermal profile with a surface density of Eg in cold 
gas and gas fraction of /g (i.e. a total mass density of E = Eg//g). Star 
formation proceeds either in a pure Kennicutt-Schmidt prescription, or the 
dynamical time variation in Eq.(Al). Cooling above 10'' K proceeds at a 
rate A and we study the simulations at several resolutions to test for conver- 
gence. 
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where numerically 
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The (vertical) dynamical time of the disk is 



^"j sech^ (^1^ mpcm ^.(21) 
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and the ratio of the dynamical time to the cooling time 

^dyn 
^cool 

L7 X 10'' ' ^ \ / Sg 



(22) 



10-22 erg crrr^ s-iy VIOM0PC-2 
The exact gravitational potential is given by 

<i>{z) = 2nGbT,J-^ log cosh 

and the pressure in hydrostatic equilibrium is 

p = nGT.gf~^bp{z) 

^g 

,lOM0pc-2 



3.3 X 10* 



0.1 



sech^ Kcm ^ . 



.(23) 

(24) 

(25) 
(26) 

(27) 



b = 



(18) 



Finally, the hydrostatic temperature for all our disks is chosen to be 
To = 10* K . (28) 

3.3 Numerical parameters and boundary conditions 

To produce simulations of a realistic ISM we make the following 
choices of parameters. In terms of resolution we must have cell 
sizes fine enough to capture the cooling of supernova rerrmants 
(Eq. 15) yet the simulation volume needs to be large enough to 
capture several scale heights of the star forming disk. In terms 
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of gas fractions and gas surface densities we cfioose values ap- 
proximating tliose in the solar neighbourhood and some varia- 
tions. In practise we chose fiducial values for the disk parameters 
(Eg = 11.61 M0pc-^ /g = 0.1) and examine this reference 
model in detail. For reference, the gas surface density of the so- 
lar neighbourhood of the Milky Way has been estimated at Eg = 
13.2 M0 pc~^, with a dynamical density of E* = 74 M© pc"^ 
(Flynn et al. 2006). 

In order to test the dependence of winds on the disk properties 
we perform a slice of the parameter space varying Eg and /g (see 
Table 1). Not all parameter combinations are explored, as we cut 
out the simulations with very small scale heights (due to resolu- 
tion constraints) and large scale heights (due to the finite box size). 
The dependence of the results on cooling, resolution, box size, star 
formation rate and run time can be seen in the Appendix. 

All our simulations were conducted in box sizes of 200 x 
200 X 1000 pc with constant cell sizes. All cells were cubic, and 
in the vertical direction the number of cells for our default reso- 
lution is 640, with corresponding cell size of 1.6 pc. We vary the 
numerical resolution using 160,320,640, 1280 cells in z, with corre- 
sponding cell sizes ranging from 6.25 — 0.78 pc. These simulations 
are denoted L2, L3, L4, L5 respectively. We also test the effect of 
adjusting our box size with simulations of 2x and 4x the width 
(see Appendix A). 

The gas surface density Eg is varied from 2.5M0 pc~^ to 
15M0 pc~^ in 8 logarithmically spaced steps followed by three 
additional steps of 30, 50 and 500 M0 pc~^ . Notably some of 
these are below the minimum surface density threshold for star for- 
mation of Schaye (2004) of 3 - 10 M0 pc"^ (although there is 
evidence that star formation proceeds below this level, e.g. Bigiel 
et al. 2008). The gas fraction /g was varied from 0.01-0.05 in 
5 logarithmic steps followed by additional steps of 0.1, 0.2, 0.5 
and 1.0. The cooling function A was varied from 3.9 x 10^'^'' to 
1.6 X 10~^^ ergcm^ s~^, and we ran additional models with the 
Sutherland & Dopita (1993) cooling function as parameterised in 
Eq. (11). Each of our experiments is evolved over 20 Myr (typi- 
cally thousands of cooling times) in order to simulate many SNe. 



4 RESULTS 

In this section we discuss the results of the simulations described in 
the previous section. We begin with a discussion of a single snap- 
shot, allowing us to investigate the instantaneous properties of the 
idealised ISM and outflow. We then move to looking at the evolu- 
tion of a simulation and the statistics we can measure before finally 
investigating the effects of all the parameters discussed in the pre- 
vious section. 

4.1 Fiducial run 

The impact of SNe depends strongly on whether they explode in 
the dense gas or in the more rarefied HIM. The supemovae in the 
disk blast bubbles in the ISM and compress the warm gas into thin 
sheets and clouds. We note that between the different simulations 
the volume of the warm medium can vary from a series of dis- 
connected, nearly spherical regions to a highly porous stratus that 
approximately covering the base of the disk potential. We will use 
the term 'clouds' to apply to both. When supernovae explode in the 
rarefied regions, either at the edge of the disk or inside previously 
evacuated bubbles, the heated gas pushes out of the central region 
and then rapidly escapes from the simulation volume in a zone of 



acceleration above and below the disk. This is the ISM portion of 
the galactic wind (i.e. the gas whose thermal energy far exceeds 
the potential barrier to escaping the disk). Some warm clouds are 
dragged along with this wind. A movie of this simulation is avail- 
able online along with time dependent versions of some of the other 
figures ^ 

In Figure 1 we show an a; — z slice of the fiducial run, at a time 
of 12 Myn We can see that the combined action of multiple SNe has 
disrupted the disk considerably, with the warm gas squeezed into 
dense sheets and globules entrained in outflowing gas, and around 
half the volume now occupied by a hot tenuous phase. The gas ap- 
pears to be in well defined phases, an HIM (greens and yellows) 
and a WNM (dark blue) with little gas at intermediate temperatures 
(sell also Fig. 2). Notably there is more temperature variation in 
the hot phase (a few orders of magnitude) than in the WNM (which 
is aU close to 10* K). The density plot also appears to show two 
distinct phases, a high and low density medivun, where the high 
densities show up in the temperature plots as WNM. In the velocity 
plot we can see a bulk vertical outflow from the disk, with veloc- 
ity correlating with height. The pressure plot shows a dramatically 
lower dynamic range than either the temperature or density plots, 
but has some distinctive shells due to individual SN remnants. The 
impression of a volume in quasi pressure equilibrium is reinforced 
by the profile plot where the temperature and density fluctuations 
appear to anti-correlate, resulting in comparatively small pressure 
variations. 

Above the plane of the disk the outflow is also very inhomo- 
geneous, containing significant turbulence as well as some warm 
clouds or globules with cometary shapes. The corresponding loca- 
tions in the density and pressure panels reveal that these clouds are 
also overdense and slightly under-pressured. In velocity the clouds 
appear to be receding from the disk at a lower velocity than the 
HIM, that is rushing past them at around 1 00 km s^ ^ . The hot wind 
is stripping the edges of these warm clouds, as evidenced by their 
tails (see also the movie online). 

After only 12 Myr the original disk has undergone consider- 
able disruption but is still observed as a connected feature in this 
slice (and the majority of the mass of the simulation remains in the 
central region). The disk has also been disrupted asymmetrically, 
with more mass pushed into the lower half space by the stochas- 
tic locations of the SNe. The externally imposed gravity will ulti- 
mately return this mass to the base of the potential, yet the com- 
bined action of the supemovae has been enough to displace it. 

Whilst we have run these simulations at different resolutions, 
it is important to note that the turbulent and chaotic nature of these 
simulations results in specific features such as individual clouds 
being at different locations or indeed absent between the different 
runs. Global properties, however, such as the outflow mass and tem- 
perature will be less stochastic, and we devote Appendix A to the 
convergence study of these properties. In general these simulations 
are numerically well converged. In the following figures we also 
include a few convergence comparisons where space allows. 

The value of the ISM pressures in our simulations are around 
10^ K cm , comparable to the pressure in simulations such as 
Joung & Mac Low (2006) and Joung et al. (2009). Estimates of 
the pressure of a star forming ISM vary, Bowyer et al. (1995) find a 
pressure of around 2 x 10** K cm~^ in the local bubble, although in 
the centre of the highly star forming region of 30 Doradus, Lopez 



^ See http: //astro . dur . ac .uk/-rmdq85 



Supernova driven winds 9 




4 5 6 7 8 -27 -26 -25 -24 -23-250-125 125 250 2 3 4 5 



loglo'r(K) login P (gcm-3) t)^ (kms-l) logiop(Kcm 3) 

Figure 1. Left to right, temperature, density, vertical velocity and pressure plots through a slice of the simulation, at time 5 Myr. Temperature is coloured 
between 10^ — 10* K, density between 20"'^^ — 10~^^ gcm~'^, from —250 to 250kms~^ and pressure from 10 — 10^ Kcm~^. On the far right 
is the profile of density, temperature and pressure along a vertical line through the centre of the slice. In dotted blue and red we show the hydrostatic density 
and pressure profiles at t = 0. Around 2; = we can see the disrupted disk in the temperature and density plots, with the warm gas squeezed into sheets and 
globules, and a significant fraction of the volume now consumed by a hot (~ 10'^"'' K) sparse phase. In the velocity plot we can see a bulk vertical outflow 
from the disk. The outflow is inhomogeneous, entraining significant turbulence as well as some warm gas, swept away from the disk. 



et al. (2011) estimate a pressure of ~ 7 x 10^ K cm~'^ from IR dust 
measurements. 

Figure 1 suggests that the hot and warm phases are quite dis- 
tinct, and we test this by inspecting the volume fractions in Fig. 2. 
The warm phase is very tightly distributed below 10*K, as we 
might expect since the only mechanism for cooling here is by adi- 
abatic expansion. The lack of intermediate temperatures suggests 
they have very short cooling times, which is consistent with a pres- 
sure equilibrium view. The hot tail of the distribution suggests the 
hottest gas either mixes with cooler gas or escapes from the simu- 
lation volume. 

Figure 3 is the density-temperature phase diagram for the fidu- 



cial model at L3 resolution (3 pc cells), which is broadly described 
by two regions. In the lower right, lying horizontally at a nearly 
constant temperature of order To = 10* K (the base of the cooling 
curve) is the WNM, which contains most of the mass. The HIM is 
in the upper left. On examination of time dependent movie of this 
simulation we see the structure in the HIM is due to multiple su- 
pernovae, each supernova blast forms a 'finger' roughly along an 
isobar, and as these shocked regions evolve and expand these lines 
descend to lower temperatures forming the mixture in the lower 
right region. As one looks to lower temperatures the fingers start 
to merge and become indistinct. We see that instantaneously we 
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Figure 2. Density and temperature probability distributions for the fiducial 
run at 10 Myr, as shown in Fig. 1, solid, dashed and dotted lines denote 
the L4, L3, L2 resolution runs, respectively. Upper panels show the mass 
fractions in temperature and density, lower panels show the corresponding 
volume fractions. We see a clear bimodality between the WNM (at low 
temperature and high density) and the HIM (at high temperature and low 
density). Almost all of the mass is in the WNM phase, but a significant 
fraction of the volume in the HIM. 
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Figure 3. Density-temperature histogram for the fiducial model at L3 
resolution. Each pixel is coloured by the fraction of cells at given p — 
T. Dashed black lines indicate lines of constant pressure, p/k-Q = 
102,10^,10* K cm ^ as indicated in the panel. We see the simulation 
volume is in an order-of-magnitude pressure equilibrium, with a bimodality 
in the gas phases into an HIM and WNM that we have segregated approx- 
imately with the dotted black line, a temperature cut at 15, 000 K. Above 
10'*K and p > 10~2*g cm"'' the cooling time of the gas is very short and 
the gas quickly cools to lO'^K. Some gas reaches lower than this tempera- 
ture due to adiabatic expansion. 



have variations in pressure within approximately one order of mag- 
nitude, and that a significant fraction of the volume is in the HIM. 



4.1.1 The characteristic temperature of the HIM 

It is interesting to consider where the characteristic temperature of 
the hot phase may appear from. We recall that the cooling function 
used in these simulations was intentionally chosen to be indepen- 



dent of temperature for T To = 10* K, and as such cannot by 
itself introduce a characteristic temperature scale, yet in Fig. 2 the 
hot gas quite clearly has a well defined peak temperature ~ lO" K. 
This is much higher than the escape temperature for the simulation 
volume (~ 10^ K, derived from Eq. 24), and as our SNe are in- 
jected just as thermal energy, there is no characteristic temperature 
for this gas. Since all of the hot gas in our simulations has been 
produced by the action of SNe it is reasonable to suppose that the 
temperature of this phase may be determined by the transition from 
the adiabatic to the momentum driven phases, as described by Cox 
(1972); Chevalier (1974) and Larson (1974). 

In this explanation, the supernovae would rapidly expand in 
the adiabatic phase until the action of cooling relative to expansion 
causes the growth of the remnant to decelerate, and the edge to form 
a cold dense shell. This shell still expands, but at a considerably 
reduced rate, driven primarily by the momentum of the shell. We 
expect the adiabatic phase to remain approximately spherical due 
to the short sound crossing time within the hot volume, however 
when the blast enters the momentum driven phase, the cooling shell 
is unstable and the remnant can become quite asymmetric. If the 
edge of the remnant reaches other sparse material the hot interior 
of the remnant can leak out (i.e. a 'chimney' such as those seen in 
Ceverino & Klypin 2009), otherwise the hot material will gradually 
be consumed into the dense shell as it radiates away its pressure 
support. 

The post shock temperature, T^, of the hot remnant at which 
the 'sag' occurs (when cooling dominates over adiabatic expan- 
sion) was calculated in Cox (1972) as 
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The obstacle which radiates away the energy of the SN is the warm 
disk gas of Fig. 1. Taking a mean density of these from Fig. 2 



3 cm 



(30) 



(p = 5 X 10~ g cm~ ) we expect a characteristic temperature 
of the remnants to be Thot ~ 3 x 10® K, very close to our HIM 
temperature of ~ 10® K. 

Another interesting application of Eq. (29) is to estimate the 
mass heated by a single supernova before it ends the adiabatic 
phase. By finding the amount of mass required to absorb the ther- 
mal energy of a supernova we derive 
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where we have neglected the initial thermal energy of the heated 
gas, the SN ejecta themselves (see also Kahn 1975), and assumed 
that none of the SN energy has yet been lost radiatively. For com- 
parison, in the model of Efstathiou (2000), a supernova evaporates 
a similar mass Afev ~ 540 Mq of cold clouds. If all this hot gas 
were to escape from the simulation without entraining any other 
material we would derive a mass loading of 
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where in Eq. (32) we have used the warm cloud density n — 
3 cm^'^ from Eq. (30), and in Eq. (33) we have used the hydro- 
static mid-plane density from Eq. (21). The mass loading is higher 
at lower surface densities (and also volume densities), at higher gas 
fractions, and for gas that cools more slowly, and increases with 
the SN energy injected. If all the gas escapes at T = Ta then this is 
an upper bound for the mass loss, since some energy will be con- 
verted to other forms such as radiation and turbulent motion, and 
for this simulation we do find the measured (3 is significantly below 
this (see section 5). Notably many versions of semi-analytic models 
such as GALFORM assume /3 close to this maximum. 

In this section we have described a snapshot of a simulation 
of a patch of the ISM with similar parameters to that of the so- 
lar neighbourhood. We have reproduced a warm and hot phase in 
order-of-magnitude pressure equilibrium, with a value similar to 
that estimated for the local volume. We have explored the relation 
between the temperature of the hot phase and related this to the den- 
sity of the warm phase via the energy of each SN and the cooling 
time of the gas. 



4.2 Time dependence 

We now turn our attention to the time dependence within our sim- 
ulation. We have seen in Fig. I that our idealised disk is disrupted 
by the energy injection from supemovae, and we are interested in 
the evolution that results from this. The injected energy can be con- 
verted into a number of forms, heating of the warm phase, the ther- 
mal energy of the hot phase, the mechanical energy of turbulence 
and the wind, the gravitational potential of the gas as it is lifted 
out of the disk, and the photons lost through radiative cooling. It 
is worth recalling that cooling is one of two ways in which energy 
can leave the simulation volume, the second being the advection of 
mass across the vertical boundaries of the simulation, taking with 
it the thermal, mechanical and gravitational potential energy of the 
gas. 

Fig. 4 is a 'space-time' plot of the onset of the outflow: time is 
along the horizontal axis, and the projected mean temperature, T, 
as a function of height is colour coded and shown on the vertical 
axis, red dots correspond to the times and location of individual SN 
injection events. In order to reduce the effects of stochastic outflows 
we performed this simulation in a larger box, of width 800 pc. The 
initially hydrostatic gas at temperature T = To seen at the far left of 
the figure is quickly replaced by gas at a range of temperatures. The 
dark blue coloured band, corresponding to T « To, episodically 
widens as a function of time, as the disk puffs up. Gas with a mean 
temperature T ~ lO'' K is seen to stream out of the disk at a range 
of velocities. From Fig. 2 we recall that there is actually very little 
gas by mass at 10® K, however by volume the mean temperature 
will be close to this. Around each supernova a plume of hot gas can 
be seen (cyan against the colder dark blue gas). At late times these 
plumes combine and drive the galactic wind. 

Comparing with the velocity lines we can see the evolution of 
the outflow velocity with time, with many structures with velocities 
in the range of 30-300 km s~ ^ . Superposed, however, are some ex- 
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Figure 4. Volume weighted mean temperature as a function of height 
and time, for the for the fiducial disk parameters yet in a wider box of 
800 X 800 X lOOOpc. At each height we have taken the average over a hori- 
zontal slice. Superimposed are red dots indicating the locations and times of 
the SN events. As the simulation progresses the activity of many SNe shock 
heat gas and drive a vertical wind from the disk at around 300kms~^. 
Dotted, dot-dashed, solid and dashed magenta lines denote outflows of 33, 
100, 300 and 1000 kms~^ respectively. Subsequent and ai'ound each su- 
pemovae can be seen a pulse (in orange) in the temperature. After a short 
(t < 1 Myr) flurry of supernova activity within the disk (2 = it53pc), the 
shocked regions begin to combine and rise out of the disk and the simulation 
volume. Occasionally, individual supernovae high above the disk (where the 
gas density is low) make a significant individual contribution to the wind. 
The 1000 km s~ ^ line has been offset to start at 6 Myr to be compared with 
the propagation of one of such temperature pulses. 



tremely steep (w.r.t. time, i.e. high velocity) discontinuities where 
much of the simulation volume rapidly experiences an increase in 
temperature. These appear to propagate from individual SNe, and 
race away from the disk with velocities in excess of 1000 km s~^, 
consistent with a sudden pressurisation of the hot phase of the ISM"^ 
This increased pressure causes stripping from the warm material as 
shocks drive in to the warmer region of the cloudy medium, adding 
to the mass of the hot phase. 

To analyse our simulations we reduced our data set down to 
the following parameters, listed below. These are chosen to give us 
a broad overview of the evolution of the star forming disk, rather 
than information on the individual cells and clouds. For these pa- 
rameters there is some freedom of definition, e.g. when one at- 
tempts to measure the pressure one could take the mid-plane pres- 
sure, the pressure within the star forming scale height b, the mean 
pressure within the simulation volume, or the mean pressure within 
a volume adjusted by some measure of the current disk scale height. 
In all cases we have attempted to choose a definition which strikes 
the balance between reducing stochasticity (some candidate mea- 
sures show considerably more noise than others) and ease of phys- 
ical intuition. 

(i) Mass ejection, Scj(t), is the amount of gas ejected from 
the disk per unit area. This is calculated from the mass advected 
through the boundary at 2; = ±500pc, divided by the surface area 
of the simulated column. This quantity is used in the calculation of 



For reference, the temperature that correspond to a given sound speed c 
is T = 7.3 X 10^ K (c/1000 km s-'^f. 
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the cosmologically important quantity /3 = Ecj/S* where we have 
identified the mass ejected from the idealised disk with the mass 
ejected from the galaxy. To achieve the nearest correspondence we 
try to maximise the volume we are measuring the loss from, i.e. the 
entire simulation volume. The corresponding normalised quantity 
is the fraction of gas remaining in the disk, /e = 1 — Eoj/Eg. 

(ii) Cold gas/Hot gas surface density is the remaining cold/hot 
gas surface densities in the simulation volume, and in combination 
with the mass ejected, sum to the initial gas surface density Eg. 

(iii) Cold volume fraction, /cow, is the volume fraction of cold 
gas, sometimes quoted in terms of the porosity 

P = -log /cold, (34) 

(Silk 2001). We distinguish between cold and hot phases at a cut-off 
of 2T(, (i.e. twice the lower limit of our cooling function). Though 
the choice of 2To may seem arbitrary, it is apparent from Fig. 2 
that the bi-modality of the warm and hot phases is quite strong, 
so the dependence of our results on the choice of temperature cut- 
off is rather low. Since the effectiveness of SNe in driving feed- 
back is highly suppressed in dense (and cold) regions, the volume 
filling factor largely determines the probability that an individual 
supernova will explode in the hot phase. The volume we study is 
z G [—250, 250] pc, as we are not interested in the hot gas far from 
the plane of the disk (where SNe do not occur). 

(iv) Pressure, p, is the mean pressure in the entire simulation 
volume. Hot material from the disk is ejected by a mean pressure 
gradient to the edge of the simulation volume, however the stochas- 
tic nature of supernova events creates a significant variation over 
small time scales and large spatial scales* and thus it is desirable to 
smooth the pressure estimate over as large a volume as possible. 

(v) Half-mass height, A1/2, is defined as the height where z £ 
[— Ai/2, A1/2] contains half the original gas mass of the disk, 

Ai/2 = min jz' : j ^ (p)^ dz > ^Egj . (35) 

At the start of the simulation this is related to the scale height by 
our choice of isothermal density profile, at A1/2 = |&log3. Large 
outflows will 'puff-up' the disk to greater scale heights, at late times 
this would become inconsistent with our star formation profile. 

(vi) Effective cooling rate, rjcff, is the total radiative cooling rate 
in the simulation volume divided by the mean SNe energy injection 
rate, 

f,, An^V 

VcB = — — ■■ ■ (36) 

l,,,£sNeioo(E*/100MQ)dA 

Conservation of energy implies that all of the energy not released as 
radiation must end up either in the wind or as gravitational poten- 
tial energy. Due to the discrete nature of time sampling with snap- 
shots (i.e. for many of the quantities such as cooling and we have 
instantaneous measurements of their time derivatives and not mea- 
surements of the integrated quantities themselves) there is some 
error on our estimate of the integrated quantities. Most susceptible 
is the estimate of the cooling rate: the tail of high-density gas seen 
in the density probability distribution function of Fig. 2, cools very 
rapidly, and our time sampling means its contribution to cooling is 
under-estimated. We will inevitably miss some cooling that would 
have occurred outside the simulation volume (although much of 

* The pressure equilibrium predicted by Spitzer (1956) holds over smaller 
spatial scales where the supersonic turbulence decays over the sound cross- 
ing time. 
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Figure 5. Generation of an outflow in the run in Fig. 4 as characterised by 
the evolution of normalised quantities described in (i)-(v) in the text. Af- 
ter a transient initial stage of ~ 5 Myr, gas starts to be ejected at a nearly 
constant rate of ^ 0.01 Mq Myr~^ pc~^. The dark blue line is the cu- 
mulative mass ejected per unit area, in units of O.2M0 pc~^. The porosity 
P = — log(/coid) of hot gas builds very quickly, green is 0.5 + 0.2P, 
implying a filling factor of the HIM of approximatelty 50%. The red line is 
the mean pressure, logjQ (p/lO^K cm"'') , disturbed from its initial value 
of 0.7 X 10^ K cm in the base of the disk by action of the SNe. Black 
line is the fraction of gas remaining in the simulation. The magenta line is 
the evolution of the scale height, Eq. (35), in terms of 0.1Aiy2 (*) /Ai/2(0)- 
The highly stochastic cyan line is rj^ff , the instantaneous cooling rate as a 
fraction of the mean SNe energy injection rate. During the first ~ 2 Myr 
the porosity in the simulation rapidly increases, after which the material be- 
gins to be ejected from the simulation in a relatively linear fashion. There 
are periods where the cooling rate increases dramatically by a factor ^ 10, 
which are closely related to SN energy injection events. Energy injection 
has not significantly puffed up the disk. 

this gas is tenuous and will have a long cooling time, little gas 
remains dense in the outflowing material). Nevertheless our high 
snapshot frequency run gives us energy conservation to ~ 1% and 
confidence that we can accurately measure the outflowing compo- 
nents from the low frequency runs (energy conservation in the sim- 
ulation itself is of course much better than this.) 

In Figure 5 we inspect these parameters for the simulation in 
Fig. 4. For the first ~ 2 Myr, the most notable feature is the rapid 
increase of porosity as the supernova blasts evacuate bubbles in the 
disk. The height of the disk remains approximately constant. As 
the simulation evolves, the remaining gas fraction declines (black 
curve) as gas leaves the simulation volume (blue curve). The mass 
lost from the simulation appears to be a nearly linear function of 
time at this stage, suggesting a constant outflow rate, which we 
investigate further in section 5. 

4.3 Comparison to a rarefaction zone 

A characteristic feature of both simulated and observed outflows 
(Steidel et al. 2010) is that the wind speed increases with height 
z above the disk, and it has been suggested that radiation driving 
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is the cause of this (Murray et al. 2005). Since radiation driving 
is not included in our modelling yet the outflow does accelerate, 
we suggest the following physical model. The combined effects of 
several supernova explosions cause the ISM pressure to increase 
substantially above the hydrostatic equilibrium value. If gravity is 
not dominant, this will lead to the higher pressure ISM expanding 
into the lower pressure regions above the disk. In the launch region 
of such an outflow, ID (plane-parallel) symmetry is a reasonable 
description of the geometry. A useful comparison is the behaviour 
of a rarefaction wave, where a homogeneous static gas is released 
into a sparse, pressure free zone, and for which the similarity solu- 
tion is 
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In such a flow, speed increases with height z and density decreases. 
This is distinct from the flow due to a single blast wave, since in 
the Sedov-Taylor phase density increases with distance from the 
blast, which is not the case for the disk outflow (Fig. 1). Notably 
this does not describe a steady wind, which would be the result of 
continuous energy injection. 

In a rarefaction wave, the acceleration is due to the pres- 
sure gradient in the outflow, and results in thermal energy be- 
ing converted to kinetic energy, and the asymptotic flow speed is 
I'max = 3co for 7 = 5/3. The outflowing gas above the disk is 
mainly warm ISM gas that is entrained by the hot SN bubbles that 
power the rarefaction wave. Figure 6 shows the behaviour of the 
simulation to be consistent with this model: velocity increases with 
height 2, but decreases with time at a given height in way predicted 
by the similarity solution. 

Notably the rarefaction is not a steady-state solution, and thus 
is not a good description of the time-averaged behaviour of the gas. 
Such behaviour should mimic the result of continuous energy injec- 
tion, where multiple overlapping SNe in the form of rarefactions or 
Sedov-Taylor blast waves (see e.g. Castor et al. 1975; Weaver et al. 
1977; McCray & Kafatos 1987) drive a large-scale wind. Our sim- 
ulations are sufficiently stochastic however that we shall leave this 
for future work. There will also be departures from a steady state 
solution as the disk consumes its gas, or in a real galaxy, has some 
gas inflow. 

4.4 Absorption features of galactic winds 

Steidel et al. (2010) proposes that the CII absorption line data is 
also well fit with velocities increasing with distance from the disk 
(in particular the lower panel of Fig. 24 of Steidel et al. 2010). The 
explanation above provides a physical mechanism for those mea- 
sured features. This is without the radiation and dust driven mech- 
anisms invoked by Murray et al. (2005); Martin (2005); Sharma 
etal. (2011). 

We pointed-out in Fig. 1 the multi-phase nature of the out- 
flow, as well as the fact that outflow speed depends on temperature. 
This is made more vivid in Fig. 7 in which we show mock 'absorp- 
tion lines' of gas selected in narrow temperature bins. These mock 
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Figure 6. Solid line shows the mean vertical velocity as a function of height 
for two times in the Eg = 2.5 Mq pc~^, /g = 0.01 simulation show- 
ing only the hot gas (where we have defined hot gas to be that above 
2 X lO** K). Red dotted line is a linear fit (a = 2.6) to the earlier snapshot 
(t = 2.5 Myr) which is then extrapolated to the later snapshot (blue dotted 
line). This shows the profile is evolving in an approximately self-similar 
fashion with the hot material accelerating away from the disk primarily due 
to its thermal energy being converted to kinetic energy. 



line profiles are simply the fraction of gas in a given temperature 
range, that is moving with a given velocity, as a function of veloc- 
ity, V. For the temperatures T < lO'^ K, the lines have their highest 
optical depths at t; ~ km s~^, and shapes which with vary lit- 
tle with temperature, T, and are almost symmetric in velocity. The 
line shapes broaden as the temperature increases, and for the hottest 
gas at r > 10^ K the line becomes asymmetric and the absorption 
centre is now ~ —100 km s~^. It is tempting to compare these to 
absorption line studies in outflows such as Martin (2005) in Nal 
and Weiner et al. (2009) in Mgll, however more work would be 
required to calculate corrections for the geometry and ionisation. 

Fig. 1 also shows colder clouds entrained inside the much hot- 
ter flow, with cometary-like tails where the cloudy medium is be- 
ing ablated by the hot gas rushing past. Absorption lines might 
arise from mass loading this hot flow either through conductive 
evaporation (see for example Boehringer & Hartquist 1987; Gnat 
et al. 2010) and/or through ablation (e.g. Hartquist et al. 1986). 
Fujita et al. (2009) investigated the warm clouds in axisymmetric 
2-dimensional simulations, where the clouds appear as Rayleigh- 
Taylor unstable cool shells and fragments that can explain the high 
velocity Na I absorption lines. We note that the metallicity of the 
gas phases is likely to be quite distinct, as the supemovae are both 
the origin of the heating and of the metals, and we intend to explore 
this in a subsequent paper. 



5 THE DEPENDENCE OF OUTFLOWS ON DISK 
PROPERTIES 

In the previous section we have discussed in detail the features of 
a simulation of a supernova-driven wind using a set of fiducial pa- 
rameters for the disk and supernova rate, the processes which drive 
it and the statistics that can be used to examine it. In this section we 
explore how the outflow properties vary and scale with the param- 
eters. We will use such scalings in the next section to integrate over 
a full galactic disk. 
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Figure 7. Normalised column density as a function of velocity, for gas with different different temperature {coloured lines). For low temperature absorbers 
(< 10'^ K) we to see a single peaked profile centred around the rest frame velocity of the disk. For higher temperatures absorbers, we see absorption at higher 
velocities relative to the disk, with velocity increasing with temperature. Only the > lO'^ K distribution appears to show any significant asymmetry. 



In Fig. 8 we plot average velocities above the simulation disk, 
for the simulations varying Eg and /g. There appears to be a strong 
trend in wind velocity, with wind speed increasing with increasing 
gas surface density, but decreasing gas fraction. There are no sim- 
ulations in the upper left as these would have a scale height larger 
than half the box size, or in the lower right as these would have a 
scale height less than 3 pc. 



5.1 Mass outflow 

Inspecting the ratio of mass outflow rate to star formation rate gives 
us an analogous property to that of Eq. (1), i.e. for a specific area 
on the disk 



/3= 5^ 



(39) 



which we use in our subsequent analysis. In theory every snapshot 
from our simulations contains an estimate of this P, as the mass 
outflow rate at a specific height, however this is rather stochastic, 
and as an alternative we calculate /? as a fit to several measurements 
of the integrated outflow 



/o* 



(40) 



which are easily obtained from each simulation snapshot. We fit the 
data samples {(ti, t/i)}"^j with the ramp function. 



fit)- 



0, t<to 
pt, t^to, 



(41) 



where the parameters to and /? are free variables. The motivation for 
choosing such a fit is that, whilst the ejection rate is nearly linear 
in most cases, there is a time (to) required for the system to reach 
a quasi steady state. This will not be a true steady state, in that the 
wind will eventually exhaust the supply of cold gas, however this 
occurs over a sufficiently long time-scale that the fit is a reasonable 
description for our simulations. 

The square error of this function can be analytically solved by 
finding linear regressions for the subsets of {(ti, yi)}"^-^ defined 
by {{ti, yj)}r=fc choosing the minimum k such that the linear 



regression t-intercept < tk.lf we define g{sk) as the t-intercept of 
the linear regression for Sk, then 

to = mm{tk ■■ g{sk) < tk,Sk = > (42) 

and P is the slope of this linear regression. 

Plots of the gas fraction remaining in the simulation volumes 
can be seen in Fig. 5 for the fiducial model, and for the set of simu- 
lations of varying Eg and /g in Fig. A9 in the Appendix where we 
also show the fits given by Eq. (41). 

In Fig. 9 we plot the mass loading /? as a function of gas sur- 
face density Eg. Each point represents a fit of /3(Eg) for the simula- 
tions varying Eg and /g. The first point to note is that our P values 
all lie below 4, and for a large range of our parameters /3 <C 1, 
i.e. our domain of parameter space switches from effective feed- 
back (more gas ejected than stars formed) to ineffective, where the 
amount of gas released is much smaller than that converted into 
stars. 

Based on jack-knife errors, our power law fit shows a signif- 
icant negative dependency, P ~ 6 (Eg/1 M© pc~^) ^ '^''^ 
implying that at high gas surface densities the feedback is less effi- 
cient. This could be due to a number of effects. Since a higher gas 
surface density will correspond to a deeper potential well, the es- 
cape velocity of the gas is higher. Secondly, the higher gaseous sur- 
face densities correspond to higher gas volume densities (Eq. (21)), 
resulting in shorter cooling times. 

Another notable dependency is that on the gas fraction. Some 
of the scatter seen in Fig. 9 actually depends systematically on the 
gas fraction, /g, with higher gas fractions showing consistently 
larger P's than the lower values. We explore this in Figure 10, 
where we have performed a simultaneous fit of P to both the gas 
surface density and the gas fraction. 



P = /3oSgY/g , 

where we find the values 

Po = 13 ± 10 
fi = 1.15±0.12 
u = 0.16±0.14, 



(43) 

(44) 
(45) 
(46) 



By construction the joint fit now no longer shows a systematic de- 
pendence on either Eg or /g. 
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Figure 8. Matrix view of simulations varying gas surface density (Sg) and 
gas fraction (/g), each panel showing a time averaged vertical velocity for 
the upper half plane of each simulation (i.e. the disk is at the base of each 
panel. Gas surface density increases from left to right, gas fraction increases 
from bottom to top. There appears to be a strong trend in wind velocity 
towards the lower right panels, i.e. a disk with low gas fraction but high gas 
surface density tends to generate a faster wind. 



Accounting for this shows a positive dependency of 
j0.i6±o.i4^ i.e. by holding the gas surface density constant but in- 
creasing the gas fraction (which reduces the gravitational strength, 
thus increasing the dynamical time and reducing the star formation 
rate) increases the mass loading. As with the dependence on gas 
surface density, we are effectively seeing a sub-linear dependence 
on star formation rate, as we decrease the star formation (increase 
the gas fraction), we see a less than proportionate drop in the out- 
flow rate. Again, the mechanism causing this should be a combina- 
tion of the processes for the Eg dependence, derived above. 

In Fig. 10 there is considerable scatter, especially at high gas 
fraction where a number of simulations have mass ejection rates 
considerably above the trend. This is most likely due to heavy dis- 
ruption of the disk out of the plane where the wind from subse- 
quent supernovae can eject it from the simulation volume. With 
such stochasticity the description of all the simulations with a sim- 
ple power law becomes inadequate. 

Our measured value for the exponents /i = 1.15 and v — 0.16 
that relate mass loading to gas surface density and gas fraction, 
P oc E"** /g , can be compared with the values from the model 
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Figure 9. The mass loading (5 (mass ejection rate vs. rate of star formation) 
as a function of gas surface density Eg. Each point represents a fit of /3 
(section 5) to a star formation simulations of varying Eg and /g. Red line 
denotes a power law fit with jack-knife errors, coloured symbols (red-blue) 
correspond to the simulations with gas fraction /g = 0.01 - 1.0 respectively. 
Vertical grey dashed line indicates the 3 Mq pc"'^ threshold for star for- 
mation from Schaye (2004). We see a significant negative dependency of 
/3 ~ (Eg/1 Mq pc~^) ""^'"^ ""^ on the gas surface density, which may 
be due to the larger gravitational potential or the higher rate of cooling (in- 
cun'ed by higher gas densities) or some combination of both. We also note 
that the scatter seems partially a function of /g, with higher gas fractions 
showing larger /3's than the lower (e.g. blue vs. green). 
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Figure 10. Joint dependence of the mass loading 13 on gas surface density. 
Eg, and gas fraction /g. Differently coloured curves correspond to simula- 
tions with different values of Egi = Eg/Mg pc~^, the thick red line is 
our best fit of the simulation points. We see a dependence of /9Eg'^^ on gas 
fraction, with a power law dependency of 0.16 it 0.15. Higher gas fractions 
for a given gas surface density imply a shallower potential well, explaining 
why the outflow efficiency increases with /g . 
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described in Section 4.1.1, which predicts scalings of/i = 8/ll = 
0.72 and v = 4/11 = 0.37. That model does not include gravity, 
and we suggest this is why the measured and predicted values dif- 
fer. To verify this we have performed a series of simulations with 
significantly higher star formation rate, described in the Appendix. 
This uses a slightly different parameterisation that is more com- 
monly used in cosmological simulations which introduces an extra 
dependence on the gas fraction, but the primary effect is an increase 
in star formation for the parameter range we study. In these runs, 
the energy injection rate is much higher, the volume filling factor 
of the hot phase much larger, and the outflow rates are correspond- 
ingly larger as well. Consequently the effect of gravity of the disk is 
much reduced. Fitting /3 oc /g to these runs yields fi = 0.82 
and v — 0.48, in much better agreement with the predictions of the 
simple model. 

It would be interesting to extend the model to account for 
the gravity of the disk, along the lines followed by Stringer et al. 
(201 1). Assume that the /3 of the hot gas in Eq.(33) is modified by 
an escape fraction /esc, which is equal to the fraction of material 
that has a temperature above the escape temperature of the simu- 
lation volume. Assuming the outflow has a range of temperatures, 
characterised by a Maxwell-Boltzmann distribution, and that only 
gas with T > Tesc escapes, the fraction is 

/esc = f /(T)dT (47) 




We have asstmied that Tesc <S Ta, i.e. the low energy tail of the 
distribution fails to escape. The net outflow will thus drop faster at 
high Eg oc Tesc, making the dependence of the mass-loading on 
Eg stronger, which is consistent with the higher fi ?» 1.15 we see 
in the lower SFR simulations. 



5.2 Radiative efficiency and energy partition in tlie ISM 

Whilst the mass loading of the galactic wind is one of the most 
cosmologically significant parameters to study, we would also like 
to evaluate the energy budgets and structure of the winds in our 
simulations. The energy injected by the SNe is absorbed into the 
gravitational binding energy, distributed into thermal and mechan- 
ical energy (both in the bulk motion of the wind and in turbulence 
throughout the simulation volume) and released as radiation (via 
cooling). 

The energy partition also enables us to evaluate a wind veloc- 
ity for the galaxy, which is commonly used to characterise feedback 
models for galaxy formation (e.g. Bower et al. 2012). The fraction 
of the energy that is incorporated into the wind, in combination with 
the mass loading, determines the overall wind speed for a galaxy. 
This is an important parameter in determining whether the wind 
can leave the galaxy and hence provide efficient quenching of star 
formation. 

By examining our simulations we can determine the fractions 
of energy that has been converted in to the different modes. In our 
fiducial simulation, we discover that a fraction of 87% was radi- 
ated, 4.5% was advected out of the computational volume as ther- 
mal energy, 5% as mechanical energy (with over half of this in 
the form of turbulent energy), 1% went into heating the simulation 



volume'', 1% went into turbulence in the simulation volume and a 
rather low 0.5% went into puffing-up the disk. The parameters here 
are averaged in a similar manner to the mass ejection rate, by taking 
the mean over snapshots after to (Eq. 42), i.e. in the quasi-steady 

regime. 

Summation of these quantities allows us to estimate rix (Eq. 
5), the fraction of power that is thermalised in to the outflow 

VT = Vthevm + r?mech , (49) 

i.e. the sum of the thermal and mechanical (bulk and turbulent) 
contributions, (the remainder going almost entirely in to cooling). 
This allows us to calculate an effective velocity Ves for the wind, 

2riT f -EsnEiooN 

where we have combined the equation for mass loading, /3 = 
Mwind/Mi,, and the thermalisation of supernova energy into the 
kinetic energy of the wind (77T), to find the specific energy in the 
wind (i.e an inversion of Eq. (5)). Notably this will be significantly 
higher than the wind velocities we see at the edge of our simulation 
volume because it includes the energy of the thermal and turbu- 
lent components. At larger distances from the galaxy, however, we 
expect this to be a more realistic estimate, as the thermal energy 
accelerates the wind and is converted in to the mechanical energy 
of the bulk flow. This is a consequence of our simulations focus- 
ing on the launch region of the galactic wind, and hence the wind 
has not yet reached its terminal velocity. Note that ram pressure 
from infalUng gas may be an important obstacle in slowing down, 
or even preventing the outflowing gas from escaping (e.g. Theuns 
et al. 2002). 

In Figure 11 we explore the dependence of the mass loading 
/3, the fraction of power in the outflow, tjt, and the effective wind 
velocity, Vea, as a function of the total surface density of the disk, 
E — Tjg/ fg. In terms of the mass loading we see a negative depen- 
dence on surface density, for comparison we have also included the 
power law fit from Eqs. (44-46). 

The fraction of power released in to the wind, r/T, appears 
to be correlated almost entirely with gas fraction /g, at high gas 
fractions much of the energy of star formation is simply radiated 
away, which is intuitive since the higher gas fractions will have 
shorter cooling times. For comparison we also show values of rjr = 
0.1 and 0.4, the former being the equivalent to the widely quoted 
10% efficiency in Larson 1974): we find star formation in disks to 
lie close to this value, except at very low gas fractions. 

The fall in outflow power in Fig. 1 1 at low surface densities 
can also be seen as a faU in the effective wind velocity. Here we 
have converted our sample values of 77T = 0.1, 0.4 into effective 
wind velocities using the power law fit for /3 in Eqs. (44-46). Each 
gas fraction appears to follow a line of approximately constant r?T, 
although there is some suggestion of a change in slope below E = 
10^ M0 pc"^. 



6 IMPACT OF OUTFLOWS ON GALAXY EVOLUTION 

In this section we apply our results from the previous section to 

the mass outflow from disk galaxies of different masses. We will 
assume a surface density profile for a galaxy and the use our fits 

^ Note that in a true steady state this fraction should be compensated by 
cooling. 
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Figure 11. Effective wind speed {upper panel), outflow efficiency (mid- 
dle panel) and mass loading (lower panel) as a function of total surface 
density S = Sg//g. Coloured lines with symbols are the simulations 
from figures (9-10), with values of the gas fraction /g as indicated. Dot- 
ted lines in the lower panel are the scalings from equations (44-46), plotted 
for /g = 0.01, 0.015, 0.2, 1.0 in the corresponding colours. Lines of con- 
stant efficiency, 77^ = 0.1 and 0.4 are shown in the middle panel (black 
dotted and dashed, respectively). Curves for the corresponding scaling of 
the effective wind speed for /g = 0.1 are shown in the upper panel. The 
outflow efficiency increases with surface density, as does the effective wind 
speed. 



for outflow efficiency as function of surface density, to deduce an 
overall feedback efficiency. 



6.1 Dependence on circular velocity from theoretical 
arguments 

In this section we take our measured dependencies of the mass 
loading parameter (which are derived for a patch of the ISM) and 
apply them to an entire disk galaxy by integrating over the surface 
of the disk. This will allow us to compare with feedback schemes 
considered in Cole et al. (2000); Bower et al. (2006) etc., which 
introduce a relation between circular velocity, mass loading, and 
effective wind speed. 

Our first step is to assume a model for a disk galaxy inside 
dark a matter halo where we follow Mo et al. (1998). The circular 
velocity of a spherical isothermal dark halo of mass A/200 is given 
by 



Vioo = lQGM2mH(z) 



(51) 



(Mo et al. 1998) where H{z) is the Hubble parameter as a function 
of redshift, z. Since the baryonic component can release energy via 
cooling, it can collapse further to become a rotationally supported 
disk. Observed bulge-less disks have a near exponential profile in 
luminous mass of the form 



E(r) = Eoexp(-r-/i?d) 



(52) 



with normalisation Eo and scale length Rd. The mass of the disk is 
thus given by 



Mi = / 27rE(r)rdr = 27rEoiiI . 



(53) 



The scale length Rd is controlled by the specific angular mo- 
mentum of the material forming the disk (e.g. Fall & Efstathiou 
1980). An exponential disk with constant rotation velocity Vd has 
angular momentum 



Jd = 47rEoV'di?d : 



(54) 



and if we parameterise in terms of the disk mass as a fraction of 
the halo mass, rud = Md/M2m, the circular velocity of the disk 
as a fraction of the halo's, Vd = Vd/V2oo and the specific angular 
momentum fraction of the disk jd/rud, we can infer the surface 
density normalisation to be 



Eo = 



2 MjVi 
n Jl 

im{z) 

ttG 



A" 



Jd 

rud 



mdWd^200 : 



(55) 



where A is the spin parameter of the isothermal halo in Eq. (51). 
Notably if we set Wd = 1 we recover the Mo et al. (1998) surface 
density equation, yet for real disks Vd > 1 as the contribution of 
baryons to the rotation velocities is not insignificant. 

We can now compute a mean mass loading j3 for such a galaxy, 
by evaluating 



/ 2Trl3t^,rdr 
J2nt-,rdr ' 



(56) 



where we will assume the surface density in star formation, E* 
follows the Kennicutt-Schmidt relation, 



AE". 



(57) 
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bulge. We have neglected the fact that there will not be any star 
formation far out in the disk if the gas surface density drops too low, 
as well as the presence of a bulge, where there may be Uttle gas and 
hence also little star formation. This will lead us to overestimate 
the wind in the tails of Fig. 12. 

To parameterise feedback in terms of the circular velocity, 
V200, we apply Eq.(55) and use our fiducial values of l3o,ii,v and 
/g, to find 



= /3o 



10 



8 10 



Figure 12. Fraction of the wind launched at each radii in the disk (Eq. 59), 
for a Kennicutt-Schmidt relation S* oc Eg , with n = 1.4, and assuming 
mass loading scales with gas surface density as /3 oc S^^, with ^ = 1.15 
(Eq. 45). Dotted line indicates the characteristic wind radius Rw/Rd for 
the galaxy, where the the local mass loading equals the net mass loading for 
the galaxy as a whole, /3 = M^ind/AJf*- 
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(63) 



(64) 



where we also assumed Ho = 71 km s ^ Mpc ^ (Freedman et al. 
2001). 

To convert to disk properties we can eliminate the spin param- 
eter with 



Rd 



XV200 f jd 
/2Q0H{z) \mdVd 



(65) 



Taking the dependence of mass loading on surface density 
found from our fits to the simulations, Eq. (43), then Eq. (56) can 
be integrated analytically. We re- write Eq. (43) in terms of the total 
surface density, E, and the gas fraction, /g, and obtain 



IM0PC-2 ' 



(58) 



giving a dependence on the gas fraction, oc /g The fraction of 
the wind laimched as function of radius is given by. 



d/w 
dir/Rd) 



2nRdl3ir)t4r) 

Afwind 



(n-/u) 



Rd 



which gives the differential rate of production of the star-formation 
driven wind, normalised by the total wind. This function is plotted 
in Fig. 12. At large radii the star formation is most effective at driv- 
ing a wind, but the net contribution to the galaxy outflow is limited 
by the low rate of star formation there. Conversely at small radii the 
wind is limited by the small area of the disk, and so it is at interme- 
diate radii where the local mass loading equals that of the galaxy as 
a whole. 

We can characterise this further by defining a wind radius Rw 

by 



/? = /?(S(i?„),/g) , 



that is, 7?w is that radius in the galaxy where the local mass loading, 
/3 = Eej/S*, equals the total mass loading of the entire galaxy, 
/3 = Mwind/M*. The wind radius for the galaxy is then given by. 



Setting jd / iTid and Vd as unity, however, yields a MW with a rather 
low circular velocity (155 km s"'^) and scale length considerably 
higher. 

The formation of the baryonic disk can increase the rotation 
velocities from V20Q both directly and indirectly. The baryons make 
their own contribution to the gravitational potential, and can also 
induce changes in the profile of the dark matter, for example due to 
adiabatic contraction (e.g. Mo et al. 2010). Even without baryons, 
there will be some adjustment to Vd due to the non-isothermal na- 
ture of halos (Navarro et al. 1997), i.e. a dependence on the concen- 
tration parameter. Here we will take Vd = 1-29 to give a circular 
speed of Vd = iiaVioo = 200 km s~^, similar to the value of the 
MW (Dehnen & Binney 1998; Flynn et al. 2006, but see also Reid 
exp[-{n- fj,)r/Rd] , (59) et al. 2009 that has the speed closer to 250 km s"^). 

Having set the circular speed, the disk scale length is implied 
by the specific angular momentum fraction in Eq. (65). For the 2.5 
kpc disk of Flynn et al. (2006) we set jd/nid ~ 0.42, i.e. the disk 
is preferentially formed of the low angular momentum baryons. A 
possible reason for the lower specific angular momentum is the de- 
layed collapse of baryons in the disk due to photo-heating (since 
disks grow in an inside-out manner, with the low angular momen- 
tum material accreted first), Navarro & Steinmetz (1997). 

Finally we should mention that the spin parameter of the MW 
may differ from 0.05, and indeed recent simulations that remove 
transient objects from halos have suggested halos have a smaller A 
(e.g. Belt et al. 2007), however we have made no account for this 
as it is outside the scope of this model. 

With these new parameters, the MW disk has a more realistic 
higher surface density, and Eq. (64) becomes 



(60) 



Rw — 



— log 
H \n — fj, 

S.ORd , 



Rd 



(61) 
(62) 



0.31 I g 



0.2 



where we have substituted in n = 1.4 for the exponent in the KS 
relation, and used the values for /x from Eq. (45). For the Milky 
Way, a disk scale length of Rd = 2.5 kpc gives a wind radius 
of J?w = 7.5 kpc, inside the solar radius but outside the galactic 



Vd 

200 kms-i 



Rd 

2.5 kpc 



rrid 
0.03 



Ho 
H{z) 



(66) 



The normalisation and scaling with Vd we find are somewhat 
below our expectations for supernova feedback. For a Milky- Way 
like halo, the star formation would remove less than one solar mass 
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of gas for every solar mass of stars formed (/3 ~ 0.31). Neverthe- 
less, halos with smaller circular velocities with the same disk ra- 
dius and disk mass fraction show increasingly effective feedback, 
(3 cx Vj~^'*, a similar scaling to energy conserving winds (e.g. 
Stringer et al. 2011). Note that the power-law dependence on Vd 
is somewhat stronger than the value of -1 found by Hopkins et al. 
(2011). Those authors also found an exponent of —0.5 for the de- 
pendence of mass loading on surface density, which is weaker than 
our exponent in Eq. (66) of /3 oc Whilst the agreement 

between these simulations is not particularly good, this is perhaps 
not surprising given that they are performed with some different 
physics, at different resolutions and using different hydrodynami- 
cal schemes. 

Despite the appeal of the above framework in supplying us 
with predictions for the mass loading in terms of redshift and the 
disk properties, there is a caveat here in our adjustment of jd/rrid 
and Vd to match the observed MW. Although we can derive this 
from observations for the MW, and the mechanism for this appears 
to be understood, it would be erroneous to suggest we have a con- 
sistent model for this, and current numerical simulations such as 
those of Scannapieco et al. (201 1) have yet to converge on the prop- 
erties of a disk for a single halo. Most concerning is that these quan- 
tities almost certainly have some implicit dependence on halo mass 
and thus there should be a corresponding adjustment to the scaling 
relation in Eq. (66). 



6.2 Dependence from observed data 

Given the approximate ingredients required to construct the for- 
malism of the previous section, it is interesting to ask whether we 
can parameterise our fit to the mass-loading, Eq. (63), with purely 
observational estimates, i.e. to compute the disk surface density 
with from observed disk properties, side-stepping the models of 
Mo et al. (1998). 

One particularly attractive method is to invert Eq. (53) to 
write the surface density in terms of the disk radius Rd and mass 
Md, where the latter can be estimated from the circular velocity 
of the disk with the TuUy-Fisher relation (Tully & Fisher 1977). 
A recent calibration of the baryonic TuUy-Fisher relations gives 
Md = 8 X 10^°MQ(V"max/200 km s"^)"* (Trachtemach et al. 
2009), application of which gives 



0.31 
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0.2 



200 km s- 



Rd 



2.5 kpc 



(67) 



which is very close to the relation in Eq. (66), including normali- 
sation and the _Rd scaling. The difference is in the exponent of Vd, 
and the dependence of Eq. (66) on md, which implicitly depends 
upon Vd as well. 

In principle it is possible to calculate the mass fraction in the 
disk from the stellar mass to halo mass function using an abun- 
dance matching approach, which would relate md to V2qo- A sin- 
gle power law, md oc M oc V200 is a good fit, although from Eq. 
(4) we see there is a dependence on the faint end slope of the stel- 
lar mass function (and at higher masses a broken power law may 
be more appropriate, e.g. Yang et al. 2003; Moster et al. 2010; Guo 
et al. 2010). Substituting this relation for md in Eq. (67) then yields 
/3 oc V^^-^Rd^ versus /3tf oc V^^'^Rld^ from Eq. (67) (taking 
= 1 . 15 for both). Finally we can try to eliminate the dependence 
on Rd, assuming Rd oc M^'^^, as inferred by Shen et al. (2003). 



This yields a scaling of oc V^ versus /Jtf oc V^ The dif- 
ference between these scalings is due to the discrepancies between 
the modelled and observed slope for the TuUy-Fisher relation and 
the uncertainty in modelling the disk mass fraction. 

Although both our scalings are strongly dependent on Vd, 
our p values were all in the range 0.01-4, so the change in feed- 
back acts more like a switch. At low disk circular velocities 
Vd ~ 140 km s"^ the feedback is high (1 < ;8 < 4) to and at the 
higher disk velocities the feedback shuts off, all over a relatively 
small range in Vd. 

To summarise, we have developed two approaches to anal- 
yse the mass loading for a galaxy based upon our estimates for 
the mass loading in our ISM patches. In Section 6.1 we take an 
analytic approximation to the properties of disk in their host ha- 
los which allows us to trace the feedback with redshift. This does, 
however, require us to make assumptions about the scaling of the 
gravitational contribution of the baryonic disks and the preferential 
accretion of low angular momentum baryons, neither of which are 
fully understood. Section 6.2 has bypassed these model concerns by 
parameterising the galaxies using the observed disk mass-velocity 
relation to directly apply the mass loadings. One price for this is the 
loss of the dependence on redshift and the cosmological evolution. 

Although these two approaches lead to different scalings, they 
do give a consistent normalisation for the feedback in the MW at 
redshift zero. In principle, one way to test this formalism is to ap- 
ply it in phenomenological models such as GALFORM, where 
such parameters as jd, ^^d and md are followed. We discuss this 
comparison further in the next section. 



6.3 Comparison to cosmological models 

We are now in a position to compare the outflow rate we measured 
in our high resolution simulations with values assumed in semi- 
analytic models such as GALFORM (Cole et al. 2000). The feed- 
back prescription for the original GALFORM was 



Vd 
^4ot 



(68) 



with values in the reference model of Vhot = 200 kms~^ and 
Ohot = 2.0. These models give a slope to the faint end of the 
galaxy luminosity function, a ~ —1.5. More recent models such 
as Bower et al. (2006) have used Qhot = 3.2 for a good match to 
the bj and K-band galaxy luminosity functions. These can be com- 
pared with our exponents from the previous paragraph, Ohot = 4.8 
and Qhot.TF ~ 2.5, which bracket the value used by Bower et al. 
(2006). For the normalisation, Cole et al. (2000) parameters yield 
0200 = 1.0 {/3 for a disk of Vd = 200 km s"^), whilst the Bower 
et al. (2006) parameters give ^200 « 17 (although this drops to 12 
using updated cosmological parameters, see Bower et al. 2012) as 
compared with our value of /320o = 0.31. The net mass loading for 
MW like galaxies obtained from our simulations is less than that 
assumed by (Cole et al. 2000) by about a factor 2, and considerably 
less than assumed by Bower et al. (2006). 

It is also interesting to consider whether the values of /3 should 
rise in starburst galaxies, where the star formation rate may be sig- 
nificantly above the normalisation of the Kennicutt-Schmidt rela- 
tion. Although our higher star formation rate simulations do show 
higher values of /5 (see Appendix A), this is only by a factor of 2, 
with l3 still falling at high gas surface densities. This suggests that 
the mechanism for galaxies to stay at high mass loadings is to re- 
main in a state with relatively low surface densities (e.g. Read et al. 
2006). 
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An alternative formulation of feedback in semi-analytics, sug- 
gested by Bower et al. (2012), is to attempt to match only the ob- 
servable portion of the stellar mass function rather than trying to 
match a slope that goes to arbitrarily faint galaxies. For example, a 
model with a constant wind speed (from the disk) ultimately pro- 
duces a faint end slope that is identical to that of the halo mass 
function. In an intermediate mass range, however, the effects of the 
gravitational potential causes material to be recycled back into the 
galaxy, producing a characteristic flat portion to the galaxy stel- 
lar mass function. By tuning the value of the wind speed, a nearly 
flat stellar mass function can be achieved over a restricted range. 
Although this mechanism cannot be extended to arbitrarily faint 
galaxies (which may be suppressed by other mechanism for exam- 
ple by re-ionization), it does provide a good fit to the observations 
with a constant /3 « 8 over this portion of the mass function. 

In contrast to some of the predictions of semi-analytic mod- 
els are the smaller estimates for the normahsation for mass loading 
found by hydrodynamic simulations. Oppenheimer et al. (2010) use 
a /3 = 2 and a «wind = 680 kms^^ to recreate the z — Q mass 
function. These simulations are at low resolution with the wind par- 
ticles partially decoupled from the surrounding gas, making them 
more comparable to semi-analytic models. Fully hydrodynamical 
simulations where the wind is coupled to the surrounding ISM are 
much harder to interpret. Resolution of these issues is beyond the 
scope of this paper, but better understanding of the differences be- 
tween semi-anal3ftic models and hydro simulations is clearly re- 
quired. 

In terms of the observed MW, Wakker et al. (2008) estimates 
the mass accretion rate to be 0.4 M© yr~^ from infalling high ve- 
locity clouds. If this is to be combined in a steady state model of a 
MW with non-negligible star formation, then j32m J; 0.4, so there 
is some tension between the observed star formation of the MW 
and the semi-analytic models that would reduce its baryon fraction, 
and our simulations lie nearer the low observed estimates. 

One option is that the semi-analytic models consistently over- 
estimate the /3200 required. In particular, there are significant de- 
generacies between /3200 and the exponent ahot. Moreover, many 
models assume that the wind scaling has a fixed energy efficiency 
(jit) and do not correctly account for the recapture of gas ejected 
from low mass galaxies (see Bower et al. 2012 for further discus- 
sion). It is entirely plausible that a careful search of parameter space 
may revel strongly mass dependent solutions much closer to those 
found here. 

On the hydrodynamical side, there are a number of physical 
processes that we neglected that may nevertheless be important. 
In terms of the gas phases we have included, the inhomogeneous 
metallicity will make an adjustment to the cooling, and larger scale 
effects such as a full 3-dimensional galactic potential along with 
shear and features such as bars and spiral arms will also play a role 
in shaping the ISM. However, it is not apparent why either of these 
effects will change the overall mass leaving the disk. In terms of 
the stellar populations we could explore the star formation distribu- 
tion in terms of the correlation with molecular clouds and also the 
clustering of stars, which may allow the explosions to strip more 
material, but this is unlikely because SNe are delayed sufficiently 
to diffuse out of their parent clouds. The large scale radiation field 
may provide an additional mechanism to accelerate the wind (Mur- 
ray et al. 2005; Hopkins et al. 2011), however in our simulations 
the thermal energy of the hot material in the disk already provides 
sufficient velocities to escape the disk. 

Potentially the largest discrepancy we have identified is the 
inconsistency of the distribution of SNe with the gas evolution. 



i.e. matching the scale height of star formation with the new scale 
height of the disk. It may even be possible to make the simulations 
completely self consistent by matching the star formation rate to the 
turbulent structure of the ISM, in a manner such as that envisaged 
by Krumholz & McKee (2005). 

Future simulations could also include the cold phase of the 
ISM by including radiative cooling below 10* K. On its own this 
would tend to reduce jS, since a cold phase removes material from 
the warm phase it would not directly increase the mass loading, 
however the physics of this brings in other processes such as self 
gravity, magnetic fields, and cosmic rays (which may be dominant 
at these scales). Magnetic fields in particular seem a candidate for 
entraining more material into the wind, although simulations such 
as Hill et al. (2012) do not find it to play a significant role. 

Overall, whilst we will include the above physical processes 
in future work, we suspect that these processes will not radically 
alter the mass-loading or significantly change the scalings we have 
found. 



7 CONCLUSIONS 

In this paper we have constructed numerically well converged sim- 
ulations of a simplified two-phase interstellar medium model, in 
which an initially isothermal and hydrostatic disk gets disrupted 
and heated by individual supernovae. By not simulating the cold 
phase of the ISM we avoided the need to introduce significantly 
more physical ingredients which require heavy algorithmic ap- 
proximations and/or fragile recipes. By restricting our simulation 
volume to only a small section of a disk, we achieve sub-parsec 
resolution, and are able to investigate the dependence of the out- 
flow on the parameters of the disk. We have included fixed gravity 
corresponding to our hydrostatic initial conditions, star formation 
that follows the Keimicutt-Schmidt relation, hydrodynamics and a 
cosmological cooling function. On scales outside the volume, the 
host disk galaxy for this toy model is reduced to the parameters of 
gas surface density, gas fraction and star formation efficiency nor- 
maUsed by the Kennicutt-Schmidt relation. 

Our simulations demonstrate the ability of supernovae to 
launch a galactic wind vertically from a disk, although we do 
not follow the subsequent evolution of the material in the halo. 
The supernovae create a turbulent ISM with very distinct hot and 
warm phases, due to the strong transition of the cooUng function at 
10* K. These phases exist in order-of-magnitude pressure equilib- 
rium, with the warm material squeezed into dense lumps, and the 
excess thermal energy of the hot material causing it to accelerate 
away from the disk. In section 4.3 we compare this to a rarefaction 
process, with the hot ISM escaping to an IGM which is compar- 
atively sparse and pressure-free. Such a model naturally leads an 
outflow with speed increasing with height above the disk but den- 
sity decreasing. 

The hot outflow entrains colder ISM gas from the disk, that 
may have relatively high metallicity. The hot gas rushes past this 
cloudy medium producing characteristic tails. Such interfaces may 
be the cites where lower ionisation lines are produced. In section 
4.4 we explore this further by calculating the normalised cross sec- 
tion of different temperature phases in our simulations, where we 
see the velocity distribution of the cooler gas is significantly be- 
neath that of the escaping material. 

In a given snapshot the precise features of our simulations vary 
greatly due to turbulence and the stochastic nature of supernovae, 
therefore we examine several global properties which are less sen- 
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sitive, such as the disk pressure, coohng rate as a fraction of the 
mean energy injection rate, disk scale height and mass ejection. 
These reveal a disk that rapidly evolves to higher porosity before 
reaching a state with an approximately constant mass ejection rate. 
This evolution of porosity is broadly reminiscent of the model by 
Silk (2001). 

We perform a range of simulation to investigate the depen- 
dence of the mass loading on gas surface density, gas fraction, and 
star formation efficiency, and fit the resulting trends with power 

laws. Our mass loadings lie in the range 0.01-4, suggesting a switch 
from a low to a high feedback regime at Vd ~ 140 kms"'^. We 
find little dependence on the normalisation of the star formation 
relation but a significant dependence on the gas fraction and sur- 
face density. The latter two can be combined to explain the bulk of 
the trends as depending on the total surface density of the disk. At 
high surface densities we find low mass loading and a high effective 
wind speed. At low surface densities the reverse is true, and there is 
an additional contribution due to an increase of the fraction of en- 
ergy radiated by cooling gas. In Section 4.1.1 we present a simple 
model where SNe blasts stall as they run into clouds swept-up by 
previous explosions that are so dense that they cool very efficiently 
predicts that mass loading depends on gas surface density and gas 
fraction as (3 — Ewind/S* oc Eg ^^''^ fg''^^ ■ These scalings are 
very close to those we find from simulations with high star forma- 
tion rate, P oc E~°'*^ /g and weaker (in terms of surface den- 
sity) than that for the pure Kennicutt relation, /3 oc Eg /g'^^. 
Our prediction for the mass loading in the solar neighbourhood is 
that each supernova results in an ejection of around 50 Mq of gas, 
or a ^ ~ 0.5, sUghtly above 0.3, our the average for the MW as a 
whole. 

The relationship between the wind velocity and thermalisa- 
tion efficiency exhibits a more complex relationship to the disk 
properties than that of the mass loading. The thermalisation effi- 
ciency appears to show a dependency on both the surface density 
and the gas fraction, and correspondingly the wind velocity does 
not show a straightforward power law implied from a constant effi- 
ciency model. For high surface densities and low gas fractions, an 
approximate 40% of the injected energy is converted into the out- 
flow's thermal, turbulent and kinetic energy components, although 
we will underestimate the cooling outside our simulation volume. 

We employ the scaling relation obtained from the simulations 
to calculate the net mass loading, P = M^ind/A^*, of an expo- 
nential disk galaxy with constant gas fraction. Using the Mo et al. 
(1998) scaling relation between disk and halo, we obtain a scaling 
with circular velocity of /3 oc V^"*'*, stronger than either energy 
or momentum-driven winds. Using the observed Tully-Fisher rela- 
tion we find a weaker dependence, $ oc V^^'^ . This compares well 
with recent semi-analytic models which assume Ohot £ [2.0, 3.2]. 

The normalisation of our net mass loading at redshift z = Q 
for a Milky- Way like galaxy is significantly lower than assumed 
in recent phenomenological models, although these models appear 
to have some degeneracy between the exponent and the normalisa- 
tion, which we will exploit in future work. Notably the mass load- 
ing only increases weakly with star formation rate but decreases 
strongly with surface density, so for starburst galaxies the feedback 
may be less efficient. Interestingly, our estimated normalisation is 
comparable with inferred values of outflow for the MW based upon 
the observed accretion and star formation. If indeed there is a higher 
mass loading, it will require supernovae to heat a larger mass of 
material to a lower temperature, or for the hot outflow to entrain a 
larger fraction of the warm ISM gas. 

The scaling we find sets the investigation of galaxy winds on a 



new footing, providing a physically motivated sub-grid description 
of winds that can be implemented in cosmological simulations and 
semi-analytic models. 
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APPENDIX A: CONVERGENCE AND PARAMETER FITS 

In this appendix we investigate the convergence properties of our 
simulations, along with the dependence upon some of the numer- 
ical and physical parameters such as the box size, simulation end 
time, the star formation rate, the cooling function and the energy of 
each supernovae. We also include evolution graphs showing the fits 
to to the mass loading which are central to this work. 

We begin by describing a set of simulations where we run an 
alternative star formation rate, which is compared to the main set 
of simulations in Section 5. In this parameterisation, the surface 
density of star formation is 

= O.lEg/tdyn , (Al) 

(more commonly used in cosmological simulations), which is ap- 
propriate for a marginally Toomre stable disk (Toomre 1964), i.e. 
the vertical dynamical time is close to the orbital time scale. Such 
prescriptions are discussed thoroughly in Schaye & Dalla Vecchia 
(2008), where they show that with self-regulating feedback this will 
recover an approximate Kennicutt-Schmidt relation of T.^^^. 

If we apply Eq.(Al) to the warm disk of our initial conditions, 
however, we generally have a much higher star formation rate due 
to the short dynamical time, which is equivalent to saying the HI 
disk is not Toomre stable. If we substitute in the tdyn from Eq. (22) 
we have a star formation rate of 

E^f Mekpc-Vr"', (A2) 

which we can see from the leading coefficient is an order of mag- 
nitude higher than Eq. (13), alhough there is some residual depen- 
dence on /g and Eg. To some extent, however, this simulates the 
conditions more relevant to a starburst galaxy. 

In Fig. Al we investigate the effect of altering the star forma- 
tion law from Eq.(13) to Eq.(Al), where the latter in general has 
much higher star formation rates due to the short vertical dynami- 
cal time of the disk. At low gas surface densities more simulations 
were possible due to the higher star formation rates. The best fit to 
the former was given in the main text, whilst the best fit to the latter 
is 

/3 ~ (20 ± 8)E-i° «2±°-°V°-*«±° °« . (A3) 
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Figure Al. Mass loading /3 as a function of Sg for the star formation laws 
in Eq. (13) {blue crosses) and Eq. (Al) {green triangles). Solid lines indicate 
the best fit for /g = 0.1 (see text for exact parameterisation). Dashed and 
dotted blue lines show the fit when varying the end time of the data used for 
the fit by ±3 Myr. 



The effect of increasing the star formation rate flattens the depen- 
dency on Eg and increases the dependency on /g, very close to the 
values predicted in Eq. (33), which is to be expected as gravity is 
much less important in these simulations (see also the discussion in 
Section 5). The relative insensitivity of /3 to the order of magnitude 
change in can be explained by the fact that the outflows are nor- 
malised by the star formation rate, so although those simulations 
have much higher outflows, the outflow per supernova deviates by 
a much smaller amount. The higher star formation rate runs can 
also be seen to have less scatter, as they are less susceptable to the 
Poisson noise of individual SN events. 

Fig. Al also illustrates the effect of the simulation end time 
on our estimate for the gas surface density dependency, by varying 
by ±3 Myr the final snapshot which is used to construct the fits for 
P (for the normal star formation rates). This shows little effect, a 
result of the outflow rates being (on average) very close to linear in 
these simulations. We perform a corresponding fit for the fiducial 
parameters in Fig. A2 where we simulate 100 Myr to check that the 
outflows we see are not a transient phenomena and continue after 
the 20 Myr of our simulations. The box width in this simulation 
was 200 pc, so we expect to see more stochasticity, and indeed we 
see fluctuations lasting many Myrs, such that the outflow estimated 
from a single 20 Myr window could show a deviation of a factor of 
a few. This is probably the main reason for the scatter in Fig. 9. 

We also test how well our simulations are converged w.r.t. the 
resolution by taking one of the high star formation rate runs and 
re-simulating it at the three resolutions L2, L3 and L4 (correspond- 
ing to a cell size of 6.3, 3.1 and 1.6 pc, respectively). In Fig. A3 
we show six different parameters, the fraction of surface density 
ejected, the mean pressure in the simulation volume, rate of cool- 
ing (as a fraction of the supernova energy injection), the porosity 
and scale heights of the disk and the wind velocity for three differ- 
ent resolutions. All of these properties appear to be well converged 
with respect to resolution, with the possible exception of the poros- 
ity and the disk scale height at the lowest resolution. With respect 
to the scale height it is notable that there is some error even at the 
initial snapshot due to the coarseness of the grid in this case. 

We explore the importance of cooling, both in broad terms 
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Figure A2. As for figures A3 {lower panel) and 4 {upper panel) but for a 
run of 100 Myr. In the lower panel the surface density ejected {dark blue 
line) has been scaled to 0.3 Mq pc~^. 



about the dependence on the magnitude of the cooling, and also 
upon our specific choice of cooling function. In Fig. A4 we look 
at the dependency of jS for the previous simulation on the magni- 
tude of the cooling function and for comparison we have included 
the Sutherland Sl Dopita (1993) cooling function for low metal- 
licity plasma, Eq. (11). The linear regression does indeed show a 
relationship albeit a weak one, with an exponent of —0.28. For the 
Sutherland & Dopita (1993) cooling function we have taken the 
magnitude of cooling to be that at the minimum, Eq. (12). The fit- 
ted j3 calculated using this figure is a factor of ~ 25% lower than 
that using our Heaviside cooling function using the same normali- 
sation. This is not quite as strong as the dependence suggested by 
Eq. (32), of -6/11 ^ -0.54. 

In Fig. A5 we make a further comparison between the Suther- 
land & Dopita (1993) cooling function and our flat cooling func- 
tion. We chose the run with the nearest normalisation (A = 1.5 x 
10~^*ergcm^ s~^) to that of the minimum of the Sutherland & 
Dopita (1993) cooling function (Eq. 12). We see a very similar 
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Figure A3. Numerical convergence of a high star formation rate run at res- 
olution of L2, L3, L4 (cell size of 6.3, 3.1 and 1.6 pc and shown in red, 
green and blue respectively). Upper left panel shows the fraction of gas that 
has left the simulation volume, middle left indicates the mean pressure in 
the simulation volume. Lower left panel shows the rate of cooling as a frac- 
tion of the mean supemovae energy injection rate, upper right shows the 
mean wind velocity, middle right shows the scale height of the disk and 
lower right shows the evolution of the porosity in the simulation. The red 
and green curves follow each-other closely, indicating good convergence. 
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Figure A4. Dependence of the mass loading parameter fi on the cool- 
ing rate. Blue crosses show the estimated f} for different values of A for 
Heaviside-shaped cooling function (A = 10~^^ ergcm"^ is the fidu- 
cial value). For comparison, the maroon plus indicates the value of /3 calcu- 
lated with the Sutherland & Dopita (1993) cooling function, at the minimum 
value of this function (Eq. (12)). 



phase distribution of the ISM, suggesting that the detailed struc- 
ture of the cooling function above 10* K does not play a large role 
in determining the features of the ISM. 

In Fig. A6 we have taken another high star formation rate sim- 
ulation and adjusted the energy associated with a single supernova. 
Here, we keep the average rate of energy injected per unit time 
fixed, but inject the energy in more (less) frequent explosions with 
less (more) energy. The variation between these simulations is sur- 
prisingly large: the behaviour of the ISM is indeed quite sensitive 
to how smooth or stochastic the injected energy is. 

In Figs. A7-A8 we investigate the dependence of the simula- 
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Figure AS. Upper panel, mass fraction of the gas in different temperature 
phases, solid, dashed lines refer to the A = 1.5 X 10"^* erg cm"^ s~^ 
Heaviside cooling function, and the Sutherland & Dopita (1993) cooling 
functions respectively. Lower panel shows the corresponding volume frac- 
tions. The fraction in each phase appears very similar, with the SD93 cool- 
ing function showing a slightly narrower temperature distribution of the hot 
phase by volume. 



tions on the size of the simulation volume. In Fig. A7 we adjust the 
vertical size of the simulation volume, i.e. whether increasing the 
volume to simulate more of the outflow adjusts the dynamics, for 
example by allowing some material to fall back to the disk. All pa- 
rameters are still computed for the original volume (±500 pc), only 
the simulation volume has been expanded. All the parameters ap- 
pear to be almost independent of this change. In Fig. A8 we adjust 
the horizontal size of simulation volume, where we have multiplied 
the box width by a factors of 2 and 4 respectively. The parameters 
here also show extremely good convergence, with the larger vol- 
umes generally showing less variation in values due to the reduced 
Poisson noise. The larger volumes also appear to show a marginal 
reduction in the evolution of the disk scale height. 

Finally, in Fig. A9 we have constructed equivalent graphs to 
that of Fig. 5, but now showing all simulations varying Eg and /g 
in Table 1. Each panel shows the time evolution of a single sim- 
ulation, showing the surface density, porosity, instantaneous cool- 
ing rate, disk height, mass ejected and pressure, along with a ramp 
function fit, Eq. (41)), to the mass ejection rates. We can see a 
strong evolution of the feedback from top left to lower right, i.e. 
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Figure A9. As Fig. 5, but for all simulations varying Eg and /g in Table 1. Additional black dashed lines have been added to show a fit to ejection {solid blue 
lines) given by Eq. (41). Sgi = E/IMq pc~^ runs from 1.5 to 11.6 from left to right panels, /g runs from 0.01 to 0.05 from lower to upper panels. Blue 
line, the amount of gas that has been ejected from the simulation, has now been scaled to units of 5 Mq pc"'^ for clarity. Green line shows 1 + P, red line, the 
mean pressure, black line, the fraction of gas remaining in the simulation, magenta line, the height of the disk and cyan line, the very stochastic instantaneous 
cooling rate as a fraction of the mean SNe energy injection rate. 
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Figure A6. Evolution of the simulations as a function of supernovae gran- 
ularity. Green line shows a run with 10^^ ergs per SN. Red line is the same 
simulation, but with 50 X the frequency of SNe, each releasing l/50th of 
the energy (2 X 10^^ ergs). Blue line has SNe at l/50th of the frequency, 
with 50 X the energy (5 X 10^2 ergs). 
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Figure A7. As for Fig. A3 but testing the effect of changing the vertical box 
size. Blue line is the outflow estimated for a simulation with the fiducial box 
height (500 pc), green line for 1 kpc. 



at high surface densities and low gas fractions the simulations de- 
velop much stronger mass ejection rates and pressures, and the disk 
is more heavily disrupted. Note that the mass ejection rate has not 
been normalised by the surface density, so much of the increase 
is due to the increased star formation in the higher surface density 
disks. In a couple of the high surface density panels the simulation 
has failed early although there are enough data points to perform a 
fit to the mass ejection. Although there is considerable stochasticity 
in the 2 parameter fits, they seem quite robust. 

In conclusion, these studies demonstrate that our simulations 
are effective at modelling a SN driven ISM and resilient to changes 
in numerical parameters. The exact nature of the cooling function 
exhibits little effect on the disk evolution, in fact the limiting fac- 
tor is largely the physical granularity of the discrete SNe and their 
locations in the disk. To this end reducing the scatter in our disk 
property dependencies could be acheived by taking a larger ensem- 



Figure A8. As Fig. A3, but for 3 different box sizes. Blue line indicates the 
run at the L3 resolution, green line with 2 X the box width and red line with 
4x the box width. 

ble of runs or alternatively by simulating larger disk areas, either of 
which increases the total number of SNe introduced. 



